plus sc10x
i.P06 and P30.PBS
check DAM signature
- P30.PBS and P30.LPS
check TDT difference/shift
load dependancies
cleaning up
processed obj
GEX.seur <- readRDS("../Spp1tdt.GEX.seur.sort1006.rds")
GEX.seur## An object of class Seurat
## 17900 features across 23829 samples within 1 assay
## Active assay: RNA (17900 features, 1246 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapcolor.FB <- ggsci::pal_igv("default")(49)[c(8,32,40,
34,26,33,28,
2,43,18)]
color.FBraw <- c("#FF6C91","lightgrey",color.FB)DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB,
pt.size = 0.1, group.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB,
pt.size = 0, group.by = "FB.info")
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age")
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age")
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc
filtering
GEX.seur <- subset(GEX.seur, subset= preAnno %in% c("Microglia") & DoubletFinder0.05 == "Singlet")
GEX.seur## An object of class Seurat
## 17900 features across 22261 samples within 1 assay
## Active assay: RNA (17900 features, 1246 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapDimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB,
pt.size = 0.1, group.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB,
pt.size = 0, group.by = "FB.info")
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age")
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age")
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc
subset2
GEX.seur <- subset(GEX.seur, subset= age %in% c("P30.PBS","P30.LPS"))
GEX.seur## An object of class Seurat
## 17900 features across 15696 samples within 1 assay
## Active assay: RNA (17900 features, 1246 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapcolor.FB2 <- color.FB[c(1:3,4:7)]DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB2)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB2,
pt.size = 0.1, group.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB2,
pt.size = 0, group.by = "FB.info")
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age")
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age")
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0.05, group.by = "sort_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0, group.by = "sort_clusters")
advanced
GEX.seur <- subset(GEX.seur, subset= percent.mt < 7.5 & nFeature_RNA > 1000 & nFeature_RNA < 3200 & nCount_RNA < 9000)
GEX.seur## An object of class Seurat
## 17900 features across 15536 samples within 1 assay
## Active assay: RNA (17900 features, 1246 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapDimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB2)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB2,
pt.size = 0.1, group.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FB2,
pt.size = 0, group.by = "FB.info")
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age")
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age")
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0.05, group.by = "sort_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0, group.by = "sort_clusters")
par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.info)
barplot(sl_stat,ylim = c(0,3800),
col = c("#FF6C91","lightgrey",color.FB),
main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(GEX.seur$FB.info), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+245,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)
subset2
re-clustering
GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1800)
# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2
head(VariableFeatures(GEX.seur), 300)## [1] "Saa3" "Cxcl10" "Ccl4" "Ccl3"
## [5] "Ccl5" "Marco" "Cxcl13" "Apoe"
## [9] "Ch25h" "Lcn2" "Ccl2" "Ccl7"
## [13] "Ccl12" "Spp1" "Rsad2" "Xist"
## [17] "Il1b" "Cd69" "Cxcl9" "Egr1"
## [21] "Acod1" "Slc13a3" "Cybb" "Iigp1"
## [25] "Il12b" "Cd74" "Ifit3" "Hist1h1b"
## [29] "Lpl" "Ifit2" "Lyz2" "Tnf"
## [33] "Cp" "Top2a" "Egr3" "Cd72"
## [37] "Ifit1" "Fth1" "Fpr1" "Ly6k"
## [41] "Aldh1a3" "Hba-a1" "Gpr84" "Ifitm3"
## [45] "Hba-a2" "Wfdc21" "Olfr889" "Mki67"
## [49] "H2-Ab1" "Hbb-bs" "Gm50255" "Lgals3"
## [53] "H2-Aa" "Wfdc17" "Gm26885" "Ccr7"
## [57] "AW112010" "Ifi27l2a" "Cxcl2" "Fn1"
## [61] "Ms4a6c" "Atf3" "Loxl2" "Ube2c"
## [65] "Tnfaip2" "Prc1" "Dennd5b" "Wdcp"
## [69] "H2-Eb1" "Cd52" "Gm26917" "Trap1a"
## [73] "Ifi211" "Pbld1" "Kcnk4" "Gdf15"
## [77] "Ifit3b" "Olfr1369-ps1" "Oasl1" "Isg15"
## [81] "Gadd45b" "Hmgb2" "Ptges" "Cxcl16"
## [85] "Egr2" "Ccr1" "Hmmr" "Hsp90ab1"
## [89] "Cenpf" "Ptgds" "Fabp5" "Ifi203"
## [93] "Hist1h2ap" "Apoc1" "Cd79a" "Ms4a4c"
## [97] "Tlr2" "Dqx1" "Ifi204" "Nusap1"
## [101] "Il10" "Hbb-bt" "Esco2" "Gm29773"
## [105] "Saa1" "E2f7" "Slfn5" "Cpd"
## [109] "Kif11" "Cdca3" "Ifi44" "Nr4a1"
## [113] "Lockd" "Bst2" "Hspa5" "Vpreb3"
## [117] "Gm10457" "Lilrb4a" "Kcnq1ot1" "Pclaf"
## [121] "Gm19409" "Ifi207" "Csf1" "Ifitm7"
## [125] "Cd63" "Tsix" "Pbk" "Ier3"
## [129] "Cd83" "Mir155hg" "Gm8251" "Sod2"
## [133] "Clec4a1" "Lhfpl4" "Ptprcap" "Rhox5"
## [137] "Cdkn1a" "Msr1" "Dnase1l3" "B930036N10Rik"
## [141] "Ifnb1" "Gm34455" "Plk1" "C4b"
## [145] "Gm15987" "Stmn1" "Srgn" "Spib"
## [149] "Ifi205" "Ifi47" "Ly6a" "Rrm2"
## [153] "H2-K1" "Hist1h1a" "Vgf" "Slc1a2"
## [157] "Hsp90aa1" "Dennd3" "A330032B11Rik" "Vnn3"
## [161] "Oasl2" "1700120C14Rik" "Hist1h2ae" "Plac8"
## [165] "1500002C15Rik" "Usp18" "Mmp12" "Ptn"
## [169] "Fcgr4" "Gm49756" "Ms4a6d" "Enah"
## [173] "Cd36" "Diras2" "Cd9" "Alas2"
## [177] "Cxcl14" "H2-Q7" "Sag" "D630023F18Rik"
## [181] "Gnaz" "Ifi30" "H2-D1" "C3"
## [185] "Tspo" "Dusp2" "Ccl9" "Ifitm2"
## [189] "Gm43409" "Ebf1" "Fcrla" "Spaar"
## [193] "Lrr1" "Mcemp1" "Ptcra" "Gbp2"
## [197] "Ifi213" "Lgals1" "Ifi209" "Vcam1"
## [201] "Syde1" "Fpr2" "Myl4" "Nes"
## [205] "Slc4a4" "Clec12a" "Mt1" "Bcl2a1d"
## [209] "mt-Nd1" "Jsrp1" "Tpx2" "Bend4"
## [213] "Arntl2" "4930414N06Rik" "Cenpe" "Ifi206"
## [217] "C3ar1" "Cd40" "Slfn1" "Parp14"
## [221] "Marcksl1" "Sox2" "Scg3" "Gm4951"
## [225] "Herpud1" "Gem" "B630019A10Rik" "Hspa1a"
## [229] "Ebf4" "Lrrc31" "Adgrg5" "Cd14"
## [233] "Neb" "Cxcr2" "Ifitm1" "Edn1"
## [237] "Ncl" "Hspa8" "Sdc4" "Nfkbia"
## [241] "Apob" "Id2" "Junb" "Pcgf2"
## [245] "Igkc" "Nfkbiz" "Irf7" "Ccna2"
## [249] "Cmpk2" "Fgl2" "Ctsc" "Lgals3bp"
## [253] "Aoah" "Chchd10" "AA467197" "Cst7"
## [257] "B230312C02Rik" "Slfn2" "Kif4" "Birc5"
## [261] "mt-Co3" "Ly6i" "S100a8" "Zfp300"
## [265] "Ifi208" "Mzb1" "P2ry12" "Tnip3"
## [269] "Ccl6" "Gpr18" "Ifit1bl1" "Oas3"
## [273] "Sox11" "Dnaja1" "3830403N18Rik" "mt-Atp6"
## [277] "Il4i1" "Gpr65" "Cpm" "Ecrg4"
## [281] "Ehf" "Gnb4" "She" "Neurod6"
## [285] "Gm44684" "Retn" "Hp" "Ppfia2"
## [289] "Olfr1393" "Mfap4" "Gm4755" "Ccdc88c"
## [293] "Dlg5" "Mov10l1" "Sst" "Cd200"
## [297] "Rit2" "Ccdc87" "Gm1141" "Col5a1"PCA
# exclude MT genes and more
# add sex-related Xist/Tsix
DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|Egr1|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AC|^AI|^AA|^AW|^AY|^BC|^Gm|^Hist|Rik$|-ps|Xist|Tsix|^Ifi|^Isg|^Mcm",
rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)
VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
c(MT_gene,
DIG,
CC_gene) )
GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))## PC_ 1
## Positive: Ccl12, Ccl2, Bst2, Slfn2, Ms4a6c, Gpr84, Lgals3bp, Ctsc, Fth1, Fcgr4
## Tspo, Srgn, Rtp4, Nme2, Calr, Fcgr1, Slfn5, Ms4a6d, Sdf2l1, Naaa
## Cxcl16, Manf, Ly6e, C5ar1, Ctsl, Cybb, Trim30a, Cd72, Oasl2, Clec2d
## Negative: P2ry12, Gpr34, Ctsd, Hpgd, Cd9, Ltc4s, Arhgap5, Sox4, Slc2a5, Pmp22
## Crybb1, Havcr2, Pmepa1, Slc40a1, Klhl24, Sgk1, Lst1, Ddit4, Il7r, Fam102b
## Thrsp, Ramp1, Slc12a2, Ttc28, Zbtb20, Tnfrsf17, Cbfa2t3, Ccl6, Crebrf, Tjp1
## PC_ 2
## Positive: Ch25h, Cybb, Ccl3, Ccl4, Aoah, Cp, Ccl5, Fpr1, Saa3, Cd72
## Cxcl10, Cd69, Ehd1, Acod1, Clec4a1, Marco, Vcam1, Ccl7, Lilrb4a, Cxcl16
## Rsad2, Oasl1, Gpr18, Gadd45b, Fpr2, Tnfaip2, Sod2, Mcemp1, Bcl2a1d, Tlr2
## Negative: Spint1, Il4ra, Gas6, Ctsl, Ccr5, Ly6e, Gbp7, Slfn2, Atp6v0a2, Gadd45g
## Naaa, Fcgr2b, Sdf2l1, Fcgr1, Cd244a, Calr, Cdc42ep3, Ccl12, Lifr, Cct6a
## Trim30a, Nav3, Tuba1a, Manf, Bhlhe41, Srgn, Nme2, Socs3, Igfbp6, Pdia3
## PC_ 3
## Positive: Rsad2, Iigp1, Usp18, Oasl2, Slfn5, Oasl1, Cxcl10, Cmpk2, Irf7, Fgl2
## Herc6, Tnfsf10, Parp14, Gbp2, Samd9l, Phf11d, Slfn9, Rnf213, Gbp5, Tlr3
## Ddx58, Samhd1, Slfn8, Stat1, Phf11b, Tgtp2, Igtp, Irgm1, Stat2, Sp100
## Negative: Npm1, Msr1, Ehd1, Eef1g, Ncl, Ran, Lpl, Anp32b, Gnl3, Tmsb4x
## Lilrb4a, Nop56, Eif5a, Cpd, Tlr2, Ftl1, C5ar1, Gpr84, Nme2, Ranbp1
## Sod2, Nop58, Pilrb1, Nme1, Tnf, Bola2, Il1b, Ctsz, Fpr1, Eif2s2
## PC_ 4
## Positive: Fcer1g, Tmsb4x, Fth1, Ly86, Ftl1, B2m, Apoe, Lyz2, Cd52, Tspo
## H2-D1, Ctsl, Cd63, Bst2, Nme2, Ssr4, Tmem176a, Actb, Prdx5, Tmem176b
## H2-K1, Psme1, Atp5g1, Zbp1, Psme2, Cst7, Spp1, Gapdh, Ly6a, Cox7b
## Negative: Macf1, Kcnq1ot1, Dst, Zeb2, Filip1l, Arhgap5, Tnf, Fam102b, Hivep3, Ccl4
## Pou2f2, Fchsd2, Ranbp2, Trip11, Slc15a3, Ankrd11, Il1b, Zfp36l1, Nav3, Tnfaip2
## Tlr2, Atp2a2, C3ar1, Baz1b, Ncl, Runx1, Tjp1, Cep350, Btg2, Sik2
## PC_ 5
## Positive: Cebpd, Il1b, Ccl12, Herpud1, Tnf, Ly6e, Ch25h, Fcgr1, Cebpb, Socs3
## C5ar1, Manf, Pdia6, Gpr84, Ms4a6d, Tmem176b, Slc2a5, Tnfaip2, Ms4a6b, Srgn
## Ms4a6c, Ccl2, Il1a, Ehd1, Sdf2l1, Saa3, Tuba1a, Gpr18, Fam102b, Nfkbia
## Negative: Apoe, Ctsb, Chd4, Ptma, B2m, Ftl1, Lpl, Gatm, Pkm, Cd52
## Stmn1, Lgals3, Itgax, Npnt, Scd2, Axl, Cstb, Mt1, Abcg1, Emp3
## Igfbp4, Lyz2, H2-D1, Atrx, Macf1, Actb, Npm1, Aplp2, Nap1l1, Lmnb1#grep("tdTomato|Spp1|Cre",VariableFeatures(GEX.seur),value = T)length(VariableFeatures(GEX.seur))## [1] 1437head(VariableFeatures(GEX.seur),500)## [1] "Saa3" "Cxcl10" "Ccl4" "Ccl3" "Ccl5"
## [6] "Marco" "Cxcl13" "Apoe" "Ch25h" "Lcn2"
## [11] "Ccl2" "Ccl7" "Ccl12" "Spp1" "Rsad2"
## [16] "Il1b" "Cd69" "Cxcl9" "Acod1" "Slc13a3"
## [21] "Cybb" "Iigp1" "Il12b" "Cd74" "Lpl"
## [26] "Lyz2" "Tnf" "Cp" "Egr3" "Cd72"
## [31] "Fth1" "Fpr1" "Ly6k" "Aldh1a3" "Gpr84"
## [36] "Wfdc21" "Olfr889" "H2-Ab1" "Lgals3" "H2-Aa"
## [41] "Wfdc17" "Ccr7" "Cxcl2" "Fn1" "Ms4a6c"
## [46] "Atf3" "Loxl2" "Tnfaip2" "Prc1" "Dennd5b"
## [51] "Wdcp" "H2-Eb1" "Cd52" "Pbld1" "Kcnk4"
## [56] "Gdf15" "Oasl1" "Gadd45b" "Ptges" "Cxcl16"
## [61] "Egr2" "Ccr1" "Ptgds" "Fabp5" "Apoc1"
## [66] "Cd79a" "Ms4a4c" "Tlr2" "Dqx1" "Il10"
## [71] "Esco2" "Saa1" "E2f7" "Slfn5" "Cpd"
## [76] "Nr4a1" "Lockd" "Bst2" "Lilrb4a" "Kcnq1ot1"
## [81] "Pclaf" "Csf1" "Cd63" "Pbk" "Ier3"
## [86] "Cd83" "Mir155hg" "Sod2" "Clec4a1" "Lhfpl4"
## [91] "Ptprcap" "Rhox5" "Cdkn1a" "Msr1" "Dnase1l3"
## [96] "Ifnb1" "Plk1" "C4b" "Stmn1" "Srgn"
## [101] "Spib" "Ly6a" "H2-K1" "Vgf" "Slc1a2"
## [106] "Dennd3" "Vnn3" "Oasl2" "Plac8" "Usp18"
## [111] "Mmp12" "Ptn" "Fcgr4" "Ms4a6d" "Enah"
## [116] "Cd36" "Diras2" "Cd9" "Alas2" "Cxcl14"
## [121] "H2-Q7" "Sag" "Gnaz" "H2-D1" "C3"
## [126] "Tspo" "Dusp2" "Ccl9" "Ebf1" "Fcrla"
## [131] "Spaar" "Lrr1" "Mcemp1" "Ptcra" "Gbp2"
## [136] "Lgals1" "Vcam1" "Syde1" "Fpr2" "Myl4"
## [141] "Nes" "Slc4a4" "Clec12a" "Mt1" "Bcl2a1d"
## [146] "Jsrp1" "Bend4" "Arntl2" "C3ar1" "Cd40"
## [151] "Slfn1" "Parp14" "Marcksl1" "Sox2" "Scg3"
## [156] "Herpud1" "Gem" "Ebf4" "Lrrc31" "Adgrg5"
## [161] "Cd14" "Neb" "Cxcr2" "Edn1" "Ncl"
## [166] "Sdc4" "Nfkbia" "Apob" "Id2" "Pcgf2"
## [171] "Nfkbiz" "Irf7" "Ccna2" "Cmpk2" "Fgl2"
## [176] "Ctsc" "Lgals3bp" "Aoah" "Chchd10" "Cst7"
## [181] "Slfn2" "Kif4" "Ly6i" "S100a8" "Zfp300"
## [186] "P2ry12" "Tnip3" "Ccl6" "Gpr18" "Oas3"
## [191] "Sox11" "Il4i1" "Gpr65" "Cpm" "Ecrg4"
## [196] "Ehf" "Gnb4" "She" "Neurod6" "Retn"
## [201] "Hp" "Ppfia2" "Olfr1393" "Mfap4" "Ccdc88c"
## [206] "Dlg5" "Mov10l1" "Sst" "Cd200" "Rit2"
## [211] "Ccdc87" "Col5a1" "Sel1l3" "Mmd2" "Plbd1"
## [216] "Spock3" "Alox8" "Mdga2" "Chrm3" "Synpr"
## [221] "Ubash3a" "Ank" "Milr1" "Ehd1" "Il2rg"
## [226] "Rnf213" "Zbtb20" "Adra2a" "Herc6" "Mrc1"
## [231] "Ctla2a" "Clec2d" "Gpr34" "Manf" "Samd9l"
## [236] "Trim17" "Prdx5" "Il1a" "Xdh" "Lox"
## [241] "Akap12" "Il1rn" "Ly86" "Tpr" "Phf11b"
## [246] "Fam20c" "Syne1" "Usp26" "Enpp2" "Cobll1"
## [251] "Ocstamp" "Robo2" "Calr" "Ctsl" "Klf2"
## [256] "Npm1" "Pcp4l1" "Ptprz1" "Smc2" "Phf11d"
## [261] "Gbp5" "Cnn3" "Spdya" "Sdf2l1" "Gbp7"
## [266] "Sox4" "Pilra" "Gbp4" "Apba2" "Ahnak2"
## [271] "Gap43" "Insl6" "Slfn9" "Olfr1423" "Rtp4"
## [276] "Igsf6" "Plcb1" "Dpp10" "Tmeff1" "A3galt2"
## [281] "Tcea3" "Kif17" "Srrm4" "Cttnbp2" "Napsa"
## [286] "Tulp2" "Zfp941" "Olfr921" "Scml4" "Gipc3"
## [291] "Cfap54" "Rflnb" "Olfr113" "Actn3" "Myof"
## [296] "Rab33a" "Gpc4" "Brs3" "Fndc3c1" "Sgk1"
## [301] "Ptger4" "Ms4a7" "Mif" "Gnao1" "C5ar1"
## [306] "Bicd1" "Anks1b" "Eif2ak2" "Stat1" "Rad50"
## [311] "Morc4" "Aox2" "Pm20d1" "Nfasc" "Zfp385b"
## [316] "Tshz2" "Bhlhe22" "Igsf3" "Synpo2" "Dkk2"
## [321] "Ptprf" "Pcdh7" "Lmntd2" "Pou2af1" "Ucn2"
## [326] "Gli1" "Slc13a5" "Ccdc92b" "Slfn4" "Pcdh8"
## [331] "Pou4f1" "Csdc2" "A4galt" "Lipo2" "Nyx"
## [336] "Lpar4" "Lst1" "Il7r" "Cfb" "Trim30a"
## [341] "Brca1" "Plp1" "Ccdc163" "Meis2" "Poln"
## [346] "H2-Q6" "Nme2" "Adgre1" "Ier2" "Naaa"
## [351] "Btg2" "Pdlim7" "Actr3b" "Spock1" "Bcl2a1a"
## [356] "Ankrd12" "Ms4a6b" "Cep290" "Rnf217" "Strip2"
## [361] "Nrgn" "Clu" "Smc1b" "Ccdc34" "Ikbke"
## [366] "Cox6a2" "Ctsb" "Taf1d" "Ftl1" "Itgal"
## [371] "Nhs" "Ctsd" "Eif5b" "Cct6b" "Ccl8"
## [376] "Crybb1" "Upk1b" "Plaur" "Icam4" "Iqch"
## [381] "Tceal3" "Pmp22" "Slc31a2" "Ccdc88a" "Lilr4b"
## [386] "Adamts1" "Arpp21" "Aspm" "Clec2i" "Mpp4"
## [391] "Mpz" "Kcnk2" "Bmp2" "Schip1" "Fam110b"
## [396] "Galnt9" "Ffar1" "Kcnc3" "Ctrl" "Cdon"
## [401] "Tmem25" "Rbms3" "Grip1" "Slc25a21" "Arg2"
## [406] "Gpr132" "Spata31d1a" "Adra1a" "Nefm" "Pcdhb4"
## [411] "Pcdhb7" "Rhod" "Plk4" "Il4ra" "Gapdh"
## [416] "Axl" "Spats2l" "Ly6e" "Ldlr" "Zbp1"
## [421] "Map7d3" "Pilrb2" "Rcan1" "Eif3a" "Socs3"
## [426] "Cxcr4" "Nrxn1" "Socs1" "Klf4" "Sgo2a"
## [431] "Brip1" "Pilrb1" "Sgo1" "Sparcl1" "Ddx58"
## [436] "Ccnb1" "Fcgr2b" "Dmrtb1" "Mns1" "Mt2"
## [441] "Sp100" "Phlda3" "Tor3a" "Smc1a" "Mndal"
## [446] "Cspg4" "Hcar2" "Knl1" "Ankrd35" "Spata7"
## [451] "Zfp811" "Agbl2" "Rragd" "Xaf1" "Ddx21"
## [456] "Etaa1os" "Runx1t1" "Ccnd1" "Fcgr1" "Fmr1nb"
## [461] "Ankrd26" "Akap9" "Eea1" "Postn" "Cadm1"
## [466] "Socs2" "Ezh2" "Samhd1" "Slamf8" "Emp3"
## [471] "Golgb1" "Icam1" "Ccdc63" "Olfr1328" "Crlf2"
## [476] "Tmem52" "Efna5" "Pclo" "Cldn7" "Nefl"
## [481] "Pcp4" "Dusp1" "Knstrn" "Anp32b" "Tmem176a"
## [486] "Ranbp1" "Morn5" "Kcnip4" "Smco3" "Adgre5"
## [491] "Tgtp1" "Rab15" "Prrt1" "Daam2" "Kcnip2"
## [496] "Cldn34b2" "Igfbp5" "Itga8" "Ly75" "Kcnh7"DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB2) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB2)
DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)
(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(desc(PC_1) ) %>% head(40)## PC_1 PC_2 PC_3 PC_4 PC_5
## Ccl12 0.14457677 -0.058118545 -0.0018318181 0.007582612 0.082396421
## Ccl2 0.11499884 0.063417638 0.0017777238 -0.036221048 0.051076874
## Bst2 0.11287944 -0.021007628 0.0277652735 0.076082452 0.010035533
## Slfn2 0.10685377 -0.071445499 0.0401899236 0.003519736 0.014894258
## Ms4a6c 0.10574457 0.038831761 -0.0014509955 0.014226352 0.052684192
## Gpr84 0.10539918 0.024039318 -0.0540462848 -0.058538073 0.063324958
## Lgals3bp 0.10417045 -0.049830290 0.0362641497 0.039165144 -0.011137749
## Ctsc 0.10095294 -0.025918417 -0.0200639846 0.007531145 0.042178753
## Fth1 0.09876637 0.050566857 -0.0450666237 0.122084274 0.018421379
## Fcgr4 0.09828840 -0.001948747 0.0240950665 0.048042450 0.021747200
## Tspo 0.09827368 -0.048911736 0.0275754640 0.092588784 -0.002869758
## Srgn 0.09716155 -0.054021370 -0.0152596736 0.030543617 0.052736727
## Rtp4 0.09018160 -0.048992222 0.0697324906 0.034217045 0.014267289
## Nme2 0.08832746 -0.053470280 -0.0531938637 0.073291733 -0.026093787
## Calr 0.08651723 -0.059477349 -0.0139172136 -0.016849079 0.040639184
## Fcgr1 0.08586060 -0.061155353 0.0478785674 0.017247831 0.074430496
## Slfn5 0.08555578 -0.006093494 0.1361551966 0.009873710 -0.026867442
## Ms4a6d 0.08544524 0.011354597 -0.0038894255 0.026874812 0.061624398
## Sdf2l1 0.08221451 -0.062916160 -0.0151158393 0.043751016 0.048960745
## Naaa 0.08077687 -0.066256678 -0.0060345557 0.023355379 0.009547661
## Cxcl16 0.08052397 0.092534600 -0.0473348881 -0.018873252 0.008203213
## Manf 0.08045464 -0.055168054 -0.0194439437 0.013993754 0.066358700
## Ly6e 0.08004240 -0.076720221 0.0467876609 0.047706314 0.074941955
## C5ar1 0.07987303 -0.052252887 -0.0549849198 -0.040784930 0.067407259
## Ctsl 0.07982985 -0.080596382 -0.0261996879 0.088109352 -0.005454459
## Cybb 0.07939603 0.167123491 -0.0126246767 -0.005137725 0.014906403
## Trim30a 0.07689381 -0.056496199 0.0823487809 -0.000542407 0.011297798
## Cd72 0.07648303 0.125824294 -0.0388548602 0.018715671 0.010916926
## Oasl2 0.07646800 0.013691289 0.1524426284 0.034955530 -0.028014516
## Clec2d 0.07581392 0.041383796 0.0422711132 -0.013785900 0.027425998
## Pdia3 0.07481310 -0.052545826 -0.0033332953 -0.023359192 0.013708430
## Ssr4 0.07423175 -0.025851064 -0.0454374940 0.071868577 -0.006559799
## Msr1 0.07401689 0.050445457 -0.0657520967 -0.047195511 0.019780185
## Cebpb 0.07222678 -0.048136119 -0.0136726779 -0.020023584 0.074376991
## Ccdc86 0.07117643 -0.033635235 -0.0006251396 -0.014705320 0.003487167
## Zbp1 0.07115050 0.005481147 0.0840959043 0.054776163 -0.014157260
## Pdia6 0.07020103 -0.047932436 -0.0154576097 -0.005955612 0.064774634
## Rnf213 0.06963596 -0.006098082 0.1007877602 -0.015361398 -0.028253033
## Usp18 0.06959384 0.022581800 0.1554107871 0.002425291 -0.011111606
## Ch25h 0.06908844 0.202951394 -0.0223728688 -0.056313752 0.074726974
## PC_6 PC_7 PC_8
## Ccl12 0.0136698459 0.0046974571 0.064177717
## Ccl2 0.0379487022 0.0632854552 0.039182200
## Bst2 0.0070299225 -0.0188225149 -0.042249980
## Slfn2 -0.0063884225 0.0095451345 0.049317330
## Ms4a6c -0.0317392773 -0.0650964964 -0.012636281
## Gpr84 0.0692366858 -0.0025107600 0.076305163
## Lgals3bp -0.0072603718 -0.0738942966 0.002982274
## Ctsc 0.0178042744 -0.0620925892 0.009439878
## Fth1 0.0251511785 -0.0631068145 -0.035943539
## Fcgr4 -0.0150690086 -0.0597129164 0.003799063
## Tspo -0.0028361490 -0.0070424635 0.011198939
## Srgn 0.0745857920 -0.0513106858 -0.024158239
## Rtp4 -0.0044163405 -0.0002538197 0.004158800
## Nme2 -0.0025394374 0.0497368308 -0.020111354
## Calr 0.0098756177 -0.0692140482 -0.072194278
## Fcgr1 0.0004209819 -0.0124601044 0.018466400
## Slfn5 0.0235324494 -0.0090961311 0.022404268
## Ms4a6d -0.0252149058 -0.0669575592 -0.041855782
## Sdf2l1 0.0126317948 -0.0564524766 -0.092176552
## Naaa -0.0297093839 -0.0199996907 0.042660378
## Cxcl16 0.0630800838 -0.0556722309 0.015232151
## Manf 0.0100804161 -0.0376948029 -0.079910284
## Ly6e 0.0327829853 -0.0632113534 -0.010175082
## C5ar1 0.0681588607 -0.0003147902 0.032370918
## Ctsl 0.0060387839 -0.0081628936 0.058153278
## Cybb -0.0623195004 -0.0477167605 -0.033990247
## Trim30a -0.0015685118 -0.0283835209 0.008363304
## Cd72 -0.0902896994 -0.0763147539 -0.021445219
## Oasl2 0.0073778915 -0.0054387363 -0.003204502
## Clec2d -0.0069487766 -0.0038654131 0.053390801
## Pdia3 -0.0028536264 -0.0852722311 -0.073132538
## Ssr4 0.0112912699 0.0080203216 -0.031416300
## Msr1 -0.0470665789 -0.0197675261 -0.036737132
## Cebpb 0.0448126464 0.0354570713 0.035654950
## Ccdc86 0.0193367799 0.0253798051 0.001065169
## Zbp1 0.0025978752 -0.0067053024 -0.023801228
## Pdia6 0.0151204901 -0.0581104673 -0.083242720
## Rnf213 0.0076861271 -0.0184715440 -0.003721215
## Usp18 0.0121419387 0.0662070870 -0.042097205
## Ch25h -0.0944706090 0.0476184573 -0.082918648(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(PC_1 ) %>% head(40)## PC_1 PC_2 PC_3 PC_4 PC_5
## P2ry12 -0.13273384 -0.005629467 0.065485216 -0.058335671 0.0151140854
## Gpr34 -0.12544252 0.015181896 0.038844037 -0.037886241 0.0028173495
## Ctsd -0.08913450 -0.002412970 0.031616554 0.017722474 0.0255570136
## Hpgd -0.07686498 0.008751071 0.022677487 -0.015074420 -0.0211753874
## Cd9 -0.06941297 0.079331428 -0.025401577 -0.009565905 -0.0391111457
## Ltc4s -0.06347304 -0.037181944 0.012272317 0.015825692 0.0042799467
## Arhgap5 -0.06343962 -0.021620906 0.042195911 -0.087940154 0.0004442022
## Sox4 -0.06032027 0.019189153 0.002317105 -0.038065717 -0.0184256103
## Slc2a5 -0.05778217 -0.017066611 0.034648291 -0.045898963 0.0554507374
## Pmp22 -0.05631634 0.010509722 0.001787813 -0.056270124 0.0295174248
## Crybb1 -0.05363518 0.003096415 0.019754398 0.031836780 -0.0555371069
## Havcr2 -0.05241728 0.025417142 0.020895802 -0.048990107 -0.0103852553
## Pmepa1 -0.05241534 -0.010455675 0.038619220 -0.038216325 0.0121668260
## Slc40a1 -0.04987126 0.006024600 0.030980473 -0.058748445 0.0051023656
## Klhl24 -0.04734085 0.027338254 0.012989055 -0.024154349 -0.0196886328
## Sgk1 -0.04631263 0.023857655 0.004370984 -0.037970720 -0.0226033893
## Lst1 -0.04405431 0.056891358 0.006394751 0.039791783 -0.0229769178
## Ddit4 -0.04263706 0.004365984 0.021918562 -0.033052296 0.0316676731
## Il7r -0.04238146 0.010493410 0.037419550 -0.030324366 0.0310134946
## Fam102b -0.04081603 -0.013585604 0.016690231 -0.084004907 0.0468995186
## Thrsp -0.04073298 0.015527898 0.019420394 -0.015585311 0.0432943542
## Ramp1 -0.04048046 0.031555175 -0.006513382 0.017861894 -0.0680257710
## Slc12a2 -0.04007155 0.017772048 0.005035451 -0.037857276 -0.0578854071
## Ttc28 -0.03850200 0.012653676 0.018642935 -0.041098095 -0.0195166438
## Zbtb20 -0.03793839 0.009000660 0.031279739 -0.042489195 0.0074691960
## Tnfrsf17 -0.03626810 0.020521373 0.010148565 -0.016869586 0.0140088553
## Cbfa2t3 -0.03488123 0.003437202 0.011728702 -0.021272566 -0.0170343673
## Ccl6 -0.03406132 0.025313415 -0.003239522 0.018979392 -0.0433961278
## Crebrf -0.03207459 0.018683094 0.024630248 -0.016528276 -0.0206936351
## Tjp1 -0.03200141 -0.004899833 0.022942055 -0.063871496 -0.0181313773
## Rapsn -0.03076310 0.023023012 0.005218324 0.010725474 -0.0160285699
## Upk1b -0.03022652 0.005547926 0.015543692 -0.037821669 0.0458513738
## Arid5b -0.02982933 0.019527703 0.030022740 -0.018714203 -0.0188086965
## Tent5c -0.02871037 0.012815758 0.008879275 -0.032698148 -0.0103643085
## S1pr1 -0.02846519 0.023880914 -0.004418825 -0.012039410 -0.0330202841
## Foxp1 -0.02803479 0.003406680 0.015871555 -0.015308127 -0.0132311851
## Csmd3 -0.02775643 -0.017580578 0.017069868 -0.061550180 0.0341682724
## Fry -0.02731603 0.005309651 0.019753401 -0.032771511 -0.0228183670
## Hacd4 -0.02722622 0.019081804 0.017959105 0.018999979 -0.0224948846
## Ppm1l -0.02714267 0.003988630 0.015116721 -0.031448084 0.0114350924
## PC_6 PC_7 PC_8
## P2ry12 -0.026156982 0.018240462 -0.0300755193
## Gpr34 0.026532943 -0.006372836 -0.0651009026
## Ctsd 0.190974686 -0.104487648 -0.1260501454
## Hpgd -0.031284412 0.049200663 -0.0083889914
## Cd9 0.151841316 -0.024208949 -0.1155900492
## Ltc4s -0.012277965 0.084374302 0.0178319841
## Arhgap5 -0.019241640 -0.059059787 -0.0760211461
## Sox4 0.011739737 0.020232577 -0.0432704868
## Slc2a5 -0.029081599 0.013571293 -0.0739490849
## Pmp22 0.110941775 -0.025663507 -0.1594533191
## Crybb1 -0.049237065 0.070428344 0.0243721919
## Havcr2 0.049704513 -0.063001520 -0.0621450743
## Pmepa1 -0.004498927 0.001010207 -0.0675075488
## Slc40a1 -0.007718383 -0.036141032 -0.0584239484
## Klhl24 -0.002514213 -0.030520337 -0.0049406889
## Sgk1 0.039579319 -0.002879153 -0.0440935219
## Lst1 -0.044907219 -0.001087295 0.0260939418
## Ddit4 -0.024443147 -0.002190533 -0.0586518156
## Il7r -0.009335944 -0.031865130 -0.0340257081
## Fam102b 0.008957638 -0.013137921 -0.0930676820
## Thrsp 0.024279140 -0.020148595 -0.0731248365
## Ramp1 -0.030205746 0.046670420 0.0623262694
## Slc12a2 0.011754490 -0.013417529 -0.0304803012
## Ttc28 0.009818731 -0.029149891 -0.0382097248
## Zbtb20 -0.001513918 -0.041651298 -0.0199665265
## Tnfrsf17 -0.001118765 0.021941105 -0.0322695726
## Cbfa2t3 -0.018356741 0.013471296 -0.0007691025
## Ccl6 0.051207150 -0.005939162 -0.0448838530
## Crebrf -0.008592314 -0.021871757 -0.0151609124
## Tjp1 -0.022404650 -0.022099382 -0.0401320487
## Rapsn 0.024991789 -0.011592548 -0.0407302242
## Upk1b 0.009742560 0.018709700 -0.0610980231
## Arid5b 0.010423645 -0.030036734 -0.0314682264
## Tent5c 0.060658674 -0.037930916 -0.0659988556
## S1pr1 0.063573154 -0.029272221 -0.0560909089
## Foxp1 -0.024303502 -0.003825656 -0.0012504121
## Csmd3 -0.009686347 0.021762739 -0.0748941670
## Fry -0.032936399 -0.006355581 0.0046956379
## Hacd4 -0.002810830 -0.007709384 0.0328864595
## Ppm1l -0.009549849 -0.001740995 -0.0240216106decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)
PCs <- 1:20
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15536
## Number of edges: 394924
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6553
## Number of communities: 15
## Elapsed time: 2 secondsRun UMAP/tSNE
GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 998)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 22:36:25 UMAP embedding parameters a = 0.9922 b = 1.112## 22:36:25 Read 15536 rows and found 20 numeric columns## 22:36:25 Using Annoy for neighbor search, n_neighbors = 20## 22:36:25 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 22:36:27 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpI5QH6X\filed2d056b0165
## 22:36:27 Searching Annoy index using 1 thread, search_k = 2000
## 22:36:31 Annoy recall = 100%
## 22:36:32 Commencing smooth kNN distance calibration using 1 thread
## 22:36:33 Initializing from normalized Laplacian + noise
## 22:36:33 Commencing optimization for 200 epochs, with 451104 positive edges
## 22:36:49 Optimization finishedDimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
GEX.seur$FB.new <- factor(as.character(GEX.seur$FB.info),
levels = levels(GEX.seur$FB.info)[c(9,8,6,7,5:3)])DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.new", split.by = "FB.new",
ncol = 4, cols = color.FB2[c(7,6,4,5,3,2,1)])
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB2)
check some markers
markers.mig <- c("Ptprc",#"Cd3d","Cd3e","Cd19",
"Cd74","Lyz2","Ccl4",
"Aif1","P2ry12","C1qa","Spp1",
"Top2a","Pcna","Mki67","Mcm6",
"Cx3cr1","Il4ra","Il13ra1","Fabp5",
"Hmox1","Ms4a7","Pf4","Tubb3",
"tdTomato-head","tdTomato-mid","tdTomato-tail","CreERT2"
)FeaturePlot(GEX.seur,
features = markers.mig,
ncol = 4)
DotPlot(GEX.seur, features = markers.mig, group.by = "FB.info",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
VlnPlot(GEX.seur, features = c("tdTomato-head","tdTomato-mid","tdTomato-tail","CreERT2"),
group.by = "FB.info", ncol = 2)
DotPlot(GEX.seur, features = c(markers.mig,"Ifit1","Ifit2"), group.by = "seurat_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
sort
GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(2,1,8,7,9,
14,
12,3,5,0,4,6,11,13,10))cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.new,
clusters=GEX.seur$sort_clusters),
main = "Cell Count",
gaps_row = c(4),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.new,
clusters=GEX.seur$sort_clusters)),
main = "Cell Ratio",
gaps_row = c(4),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
markers
# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "sort_clusters"
#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE,
# min.pct = 0.05,
# test.use = "MAST",
# logfc.threshold = 0.25)
GEX.markers.pre <- read.table("sc10x_LYN.marker.subset2_sort1007.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)## # A tibble: 120 x 7
## # Groups: cluster [15]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 0 0.887 0.989 0.922 0 2 Gpr34
## 2 0 0.858 1 0.973 0 2 P2ry12
## 3 5.16e-226 0.604 0.981 0.963 9.24e-222 2 Rgs10
## 4 1.17e-204 0.652 0.921 0.855 2.09e-200 2 Rnase4
## 5 4.07e-203 0.707 0.923 0.804 7.29e-199 2 Cd9
## 6 1.13e-192 0.724 0.796 0.665 2.03e-188 2 Hpgd
## 7 1.56e-143 0.603 0.586 0.48 2.80e-139 2 Ang
## 8 4.46e-130 0.613 0.583 0.477 7.98e-126 2 Lst1
## 9 4.83e-170 0.923 0.962 0.847 8.64e-166 7 H2-D1
## 10 5.05e-108 0.582 0.996 0.926 9.04e-104 7 Gpr34
## # ... with 110 more rowsGEX.markers.sort <- GEX.markers.pre
GEX.markers.sort$cluster <- factor(as.character(GEX.markers.sort$cluster),
levels = levels(GEX.seur$sort_clusters))
markers.pre_t48 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.1) %>%
top_n(n = 48, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t60 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t120 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
top_n(n = 120, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$geneDotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[1:64])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[65:128])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[129:192])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[193:256])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[257:320])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[321:384])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[385:448])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[449:498])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[499:529])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
GEX.seur$cnt <- as.character(GEX.seur$FB.info)
GEX.seur$cnt[GEX.seur$cnt %in% c("P30.LPS.TDT1","P30.LPS.TDT2")] <- "P30.LPS.TDT"
GEX.seur$cnt <- factor(GEX.seur$cnt,
levels = c("P30.PBS.CTR","P30.PBS.MIG","P30.PBS.TDT",
"P30.LPS.CTR","P30.LPS.MIG","P30.LPS.TDT"))
color.cnt <- rev(color.FB2)[c(5:7,1:2,4)]subset2 DEGs
DEGs.a
GEX.comp <- GEX.seur
Idents(GEX.comp) <- "cnt"GEX.comp## An object of class Seurat
## 17900 features across 15536 samples within 1 assay
## Active assay: RNA (17900 features, 1437 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapGEX.comp$cnt[1:6]## AAACCCAAGGTCGACA-1 AAACCCAAGGTTCTTG-1 AAACCCACAAACTCGT-1 AAACCCACACGCGGTT-1
## P30.LPS.CTR P30.LPS.TDT P30.PBS.MIG P30.LPS.TDT
## AAACCCACATCCGTGG-1 AAACCCAGTATCGAAA-1
## P30.PBS.TDT P30.LPS.MIG
## 6 Levels: P30.PBS.CTR P30.PBS.MIG P30.PBS.TDT P30.LPS.CTR ... P30.LPS.TDTDEGs.a <- list()
for(mm in c("P30.PBS","P30.LPS")){
for(nn in c("CTR","MIG")){
DEGs.a[[paste0(mm,".TDTvs",nn)]] <- FindAllMarkers(subset(GEX.comp, subset = cnt %in% c(paste0(mm,".",nn),paste0(mm,".TDT")) & age %in% c(mm)),
assay = "RNA",
only.pos = T,
min.pct = 0.05,
logfc.threshold = 0.1,
test.use = "MAST")
}
}## Calculating cluster P30.PBS.CTR##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.PBS.TDT##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.PBS.MIG##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.PBS.TDT##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.LPS.CTR##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.LPS.TDT##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.LPS.MIG##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.LPS.TDT##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!#DEGs.anames(DEGs.a)## [1] "P30.PBS.TDTvsCTR" "P30.PBS.TDTvsMIG" "P30.LPS.TDTvsCTR" "P30.LPS.TDTvsMIG"# save DEGs' table
df.DEGs.a <- data.frame()
for(nn in names(DEGs.a)){
df.DEGs.a <- rbind(df.DEGs.a, data.frame(DEGs.a[[nn]], contrast = nn))
}
#write.csv(df.DEGs.a, "DEGs.a.subset2.csv")DEG stat
cut.padj = 0.05
cut.log2FC = 0.3
cut.pct1 = 0.05
df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()## contrast cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR 8
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT 29
## 3 P30.LPS.TDTvsMIG P30.LPS.TDT 51
## 4 P30.LPS.TDTvsMIG P30.LPS.MIG 17
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR 2
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT 7
## 7 P30.PBS.TDTvsMIG P30.PBS.TDT 1
## 8 P30.PBS.TDTvsMIG P30.PBS.MIG 3pp.stat.DEG <- list()df.DEGs.a$cluster <- factor(as.character(df.DEGs.a$cluster),
levels = levels(GEX.seur$cnt)[c(4:6,1:3)])pp.stat.DEG[["a"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["a"]]
cut.padj = 0.05
cut.log2FC = 0.25#0.3
cut.pct1 = 0.05
df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()## contrast cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR 13
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT 37
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG 20
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT 72
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR 4
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT 11
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG 4
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT 3pp.stat.DEG[["a1"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["a1"]]
DEG plot
pp.DEGs.a <- lapply(DEGs.a, function(x){NA})names(DEGs.a)## [1] "P30.PBS.TDTvsCTR" "P30.PBS.TDTvsMIG" "P30.LPS.TDTvsCTR" "P30.LPS.TDTvsMIG"DEGs.a.combine <- DEGs.a
DEGs.a.combine <- lapply(DEGs.a.combine, function(x){
for(zz in c("P30.LPS.MIG","P30.LPS.CTR","P30.PBS.MIG","P30.PBS.CTR")){
x[x$cluster==zz,"avg_log2FC"] <- -x[x$cluster==zz,"avg_log2FC"]
}
x
})pp.DEGs.a$P30.LPS.TDTvsCTR <- DEGs.a.combine$P30.LPS.TDTvsCTR %>%
mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTR"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.LPS, TDT vs CTR") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,2))
pp.DEGs.a$P30.LPS.TDTvsCTR
pp.DEGs.a$P30.LPS.TDTvsMIG <- DEGs.a.combine$P30.LPS.TDTvsMIG %>%
mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("MIG"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.LPS, TDT vs MIG") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,1))
pp.DEGs.a$P30.LPS.TDTvsMIG
pp.DEGs.a$P30.PBS.TDTvsCTR <- DEGs.a.combine$P30.PBS.TDTvsCTR %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTR"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.PBS, TDT vs CTR") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,2))
pp.DEGs.a$P30.PBS.TDTvsCTR
pp.DEGs.a$P30.PBS.TDTvsMIG <- DEGs.a.combine$P30.PBS.TDTvsMIG %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("MIG"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.PBS, TDT vs MIG") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,2))
pp.DEGs.a$P30.PBS.TDTvsMIG
DEGs.b
GEX.comp <- GEX.seur
Idents(GEX.comp) <- "age"GEX.comp## An object of class Seurat
## 17900 features across 15536 samples within 1 assay
## Active assay: RNA (17900 features, 1437 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapGEX.comp$age[1:6]## AAACCCAAGGTCGACA-1 AAACCCAAGGTTCTTG-1 AAACCCACAAACTCGT-1 AAACCCACACGCGGTT-1
## P30.LPS P30.LPS P30.PBS P30.LPS
## AAACCCACATCCGTGG-1 AAACCCAGTATCGAAA-1
## P30.PBS P30.LPS
## Levels: P06 P30.PBS P30.LPSlist.b <- list(P30.CTR = c("P30.PBS.CTR","P30.LPS.CTR"),
P30.MIG = c("P30.PBS.MIG","P30.LPS.MIG"),
P30.TDT = c("P30.PBS.TDT","P30.LPS.TDT"))
list.b## $P30.CTR
## [1] "P30.PBS.CTR" "P30.LPS.CTR"
##
## $P30.MIG
## [1] "P30.PBS.MIG" "P30.LPS.MIG"
##
## $P30.TDT
## [1] "P30.PBS.TDT" "P30.LPS.TDT"names(list.b)## [1] "P30.CTR" "P30.MIG" "P30.TDT"DEGs.b <- list()
for(nn in names(list.b)){
DEGs.b[[nn]] <- FindAllMarkers(subset(GEX.comp, subset = age %in% c("P30.PBS","P30.LPS") & cnt %in% list.b[[nn]]),
assay = "RNA",
only.pos = T,
min.pct = 0.05,
logfc.threshold = 0.1,
test.use = "MAST")
}## Calculating cluster P30.PBS##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.LPS##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.PBS##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.LPS##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.PBS##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Calculating cluster P30.LPS##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!#DEGs.bnames(DEGs.b)## [1] "P30.CTR" "P30.MIG" "P30.TDT"# save DEGs' table
df.DEGs.b <- data.frame()
for(nn in names(DEGs.b)){
df.DEGs.b <- rbind(df.DEGs.b, data.frame(DEGs.b[[nn]], condition = nn))
}
#write.csv(df.DEGs.b, "DEGs.b.subset2.csv")DEG stat
cut.padj = 0.05
cut.log2FC = 0.3
cut.pct1 = 0.05
df.DEGs.b %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(condition,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()## condition cluster up.DEGs
## 1 P30.CTR P30.PBS 425
## 2 P30.CTR P30.LPS 512
## 3 P30.MIG P30.PBS 321
## 4 P30.MIG P30.LPS 391
## 5 P30.TDT P30.PBS 388
## 6 P30.TDT P30.LPS 423pp.stat.DEG <- list()df.DEGs.b$cluster <- factor(as.character(df.DEGs.b$cluster),
levels = c("P30.PBS","P30.LPS"))pp.stat.DEG[["b"]] <- df.DEGs.b %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(condition,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=condition, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["b"]]
cut.padj = 0.05
cut.log2FC = 0.3
cut.pct1 = 0.05
df.DEGs.b %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(condition,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()## condition cluster up.DEGs
## 1 P30.CTR P30.PBS 425
## 2 P30.CTR P30.LPS 512
## 3 P30.MIG P30.PBS 321
## 4 P30.MIG P30.LPS 391
## 5 P30.TDT P30.PBS 388
## 6 P30.TDT P30.LPS 423pp.stat.DEG[["b"]] <- df.DEGs.b %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(condition,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=condition, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["b"]]
DEG plot
pp.DEGs.b <- lapply(DEGs.b, function(x){NA})names(DEGs.b)## [1] "P30.CTR" "P30.MIG" "P30.TDT"DEGs.b.combine <- DEGs.b
DEGs.b.combine <- lapply(DEGs.b.combine, function(x){
x[x$cluster=="P30.PBS","avg_log2FC"] <- -x[x$cluster=="P30.PBS","avg_log2FC"]
x
})CTR
pp.DEGs.b$P30.CTR <- DEGs.b.combine$P30.CTR %>%
mutate(label=ifelse(((p_val_adj < 1e-100 & avg_log2FC<0) | (p_val_adj < 1e-100 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"P30.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"P30.LPS","None")),
padj=ifelse(p_val_adj<1e-300,p_val_adj+1e-300,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("P30.PBS"="Skyblue",
"P30.LPS"="pink",
"None"="grey")) +
theme_classic() + labs(title="CTR, P30.PBS vs P30.LPS") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1,2))
pp.DEGs.b$P30.CTR
MIG
pp.DEGs.b$P30.MIG <- DEGs.b.combine$P30.MIG %>%
mutate(label=ifelse(((p_val_adj < 1e-80 & avg_log2FC<0) | (p_val_adj < 1e-80 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"P30.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"P30.LPS","None")),
padj=ifelse(p_val_adj<1e-300,p_val_adj+1e-300,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("P30.PBS"="Skyblue",
"P30.LPS"="pink",
"None"="grey")) +
theme_classic() + labs(title="MIG, P30.PBS vs P30.LPS") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1,2))
pp.DEGs.b$P30.MIG
TDT
pp.DEGs.b$P30.TDT <- DEGs.b.combine$P30.TDT %>%
mutate(label=ifelse(((p_val_adj < 1e-80 & avg_log2FC<0) | (p_val_adj < 1e-80 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"P30.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"P30.LPS","None")),
padj=ifelse(p_val_adj<1e-300,p_val_adj+1e-300,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("P30.PBS"="Skyblue",
"P30.LPS"="pink",
"None"="grey")) +
theme_classic() + labs(title="TDT, P30.PBS vs P30.LPS") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1,2))
pp.DEGs.b$P30.TDT
signature score
#### 10x data, calculate signature score
#
## The code below is from Adam Hamber
## 2D scoring by Itay
get_controls <- function(counts, gene.list, verbose=F, control.genes.per.gene=10)
{
# Itay: "Such scores are inevitably correlated with cell complexity so to avoid
# that I subtract a "control" score which is generated by averaging over a control
# gene set. Control gene sets are chosen to contain 100 times more genes than the
# real gene set (analogous to averaging over 100 control sets of similar size) and
# to have the same distribution of population/bulk - based expression levels as the
# real gene set, such that they are expected to have the same number of "zeros" and
# to eliminate the correlation with complexity."
# ---------------------------------------------------------------------------------
# Going to find control points by finding the closest genes in terms of expression level and % of the time we observe it
if(verbose){cat(sprintf("Finding %s background genes based on similarity to given gene set [%s genes] \n",
control.genes.per.gene*length(gene.list), length(gene.list)))}
cat("Summarizing data \n")
summary = data.frame(gene=row.names(counts), mean.expr = Matrix::rowMeans(counts), fract.zero = Matrix::rowMeans(counts==0), stringsAsFactors = F)
#summary = data.frame(gene=row.names(counts), mean.expr = apply(counts,1,mean), fract.zero = apply(counts==0,1,mean), stringsAsFactors = F)
summary$mean.expr.s = scale(summary$mean.expr)
summary$fract.zero.s = scale(summary$fract.zero)
actual.genes = summary[summary$gene %in% gene.list,]
background.genes = summary[!summary$gene %in% gene.list,]
#find the 10 closest genes to each cell cycle marker gene and add them to the lists of control genes
get_closest_genes <- function(i)
{
background.genes$dist = sqrt((background.genes$mean.expr.s - actual.genes$mean.expr.s[i])^2 +
(background.genes$fract.zero.s - actual.genes$fract.zero.s[i])^2)
ordered = background.genes$gene[order(background.genes$dist)]
ordered = ordered[!ordered %in% controls] # don't take genes that already appear in the list
closest = head(ordered, n=control.genes.per.gene)
return(closest)
}
controls = c();
for (i in 1:length(gene.list)){
#info(sprintf("Finding %s control genes for %s", control.genes.per.gene, gene.list[i]))
closest = get_closest_genes(i)
#info(sprintf("Found %s: ", length(closest)))
controls = unique(c(controls, closest))
}
if(verbose){cat(sprintf("Control gene selection complete. %s genes found. \n", length(controls)))}
#print(controls)
return(controls)
}
## Define calculate function
calculate_signature_score <- function(count_matrix, gene_list){
control_gene <- get_controls(counts = count_matrix,
gene.list = gene_list)
signature_score <- colMeans(count_matrix[gene_list, ], na.rm = TRUE) -
colMeans(count_matrix[control_gene, ], na.rm = TRUE)
return(signature_score)
}
add_geneset_score <- function(obj, geneset, setname){
score <- calculate_signature_score(as.data.frame(obj@assays[['RNA']]@data),
geneset)
obj <- AddMetaData(obj,
score,
setname)
return(obj)
}## proc_DEG()
# to process edgeR result for DEGs' comparison
# mat_cut is a matrix after advanced filtering, 'like TPM > n in at least one condition'
#
proc_DEG <- function(deg, p.cut=0.05, FC.cut = 2, padj=TRUE, abs=TRUE, mat_cut=NULL, gene_cut=NULL){
rownames(deg) <- deg$gene
if(padj==TRUE){
deg <- deg %>% filter(padj < p.cut)
}else{
deg <- deg %>% filter(pvalue < p.cut)
}
if(abs==TRUE){
deg <- deg %>% filter(abs(FC) > FC.cut)
}else if(FC.cut >0){
deg <- deg %>% filter(FC > FC.cut)
}else{
deg <- deg %>% filter(FC < FC.cut)
}
if(!is.null(mat_cut)){
deg <- deg[rownames(deg) %in% rownames(mat_cut),]
}
if(!is.null(gene_cut)){
deg <- deg[rownames(deg) %in% gene_cut,]
}
return(deg)
}#saveRDS(GEX.seur, "Spp1tdt.subset2_P30LPS_P30PBS.sort1007.rds")DEG.list <- list(LPS.TDT= (df.DEGs.a %>% filter(p_val_adj < 0.05 &
abs(avg_log2FC) > 0.1 &
pct.1 > 0.05 & contrast == "P30.LPS.TDTvsCTR" &
cluster == "P30.LPS.TDT"))$gene,
LPS.CTR= (df.DEGs.a %>% filter(p_val_adj < 0.05 &
abs(avg_log2FC) > 0.1 &
pct.1 > 0.05 & contrast == "P30.LPS.TDTvsCTR" &
cluster == "P30.LPS.CTR"))$gene)DEG.list## $LPS.TDT
## [1] "Ly86" "Xist" "Cd52" "Fn1"
## [5] "Apoe" "tdTomato-mid" "H2-D1" "H2-K1"
## [9] "Cd72" "Gatm" "Fth1" "Ctss"
## [13] "Gm21860" "Ms4a6d" "B2m" "Lyz2"
## [17] "H2-Q7" "E230029C05Rik" "Cybb" "Cd74"
## [21] "Ms4a6c" "Trps1" "Cp" "Spp1"
## [25] "Ifi204" "C4b" "Rps12" "Fcgr4"
## [29] "Rps29" "Ccl5" "Ifitm3" "Ifi211"
## [33] "Mpeg1" "Cd180" "Saa3" "Tspo"
## [37] "Milr1" "Ifi207" "Sash1" "Chst1"
## [41] "Ch25h" "mt-Nd1" "H2-Q6" "Ifi27l2a"
## [45] "Axl" "Npc2" "Clec4a1" "AW112010"
## [49] "Ctsz" "Lgals3bp" "Ccl2" "Nceh1"
## [53] "Rpl39" "Rplp1" "Ifi213" "mt-Atp6"
## [57] "Ms4a4c" "Xdh" "Malat1" "Itm2b"
## [61] "Lilrb4a" "Lst1" "Slfn5" "Aoah"
## [65] "Rps20" "Bst2" "Ifitm2" "Gpr65"
## [69] "mt-Co2" "Zbp1" "P2rx4" "Fyb"
## [73] "Epb41l3" "Ptger4" "Eif2s3y" "Clec2d"
## [77] "H2-Ab1" "Fcer1g" "Rps28" "Lair1"
## [81] "Rps27" "Clec4a3" "Cpd" "mt-Cytb"
## [85] "Rps24" "Ifi209" "Slc6a6" "Rps3"
## [89] "Rps21" "Hcar2" "Tor3a" "Cd48"
## [93] "Itgam" "Iigp1" "Tmcc3" "Rpl37a"
## [97] "Tlr2" "Ly9" "Ifi206" "Zeb2"
## [101] "Srgap2" "AU020206" "mt-Co3" "H2-M3"
## [105] "Ifit3" "Ptprc" "Unc93b1" "Fabp5"
## [109] "Rsad2" "Cd14" "Ccl9" "Prdx5"
## [113] "Gm4951" "Psap" "Ctsc" "Mt1"
## [117] "Tnfsf13b" "Sp100" "Rpl23" "Rasa4"
## [121] "Anxa5" "Ccl3" "Tmem86a" "Rpl10a"
## [125] "Ccl6" "Plgrkt" "Slc31a2" "Mbnl1"
## [129] "Fgl2" "Oasl2"
##
## $LPS.CTR
## [1] "tdTomato-head" "Hspa8" "Dnaja1" "Hsp90ab1"
## [5] "Hsp90aa1" "Cd164" "Fcrls" "Tmem119"
## [9] "Hspa1a" "Hspa1b" "Hspa5" "Hsph1"
## [13] "Cacybp" "Hspd1" "Hspe1" "G530011O06Rik"
## [17] "Pmepa1" "Fkbp4" "Lrba" "Sparc"
## [21] "Tuba1a" "Basp1" "Stip1" "Mat2a"
## [25] "Gt(ROSA)26Sor" "Abi3" "Chordc1" "Slc2a5"
## [29] "Manf" "Ubb" "Cct3" "Cebpd"
## [33] "Eif1" "Cd34" "P4ha1" "Ptges3"
## [37] "Ubc" "Cct7" "Cct6a" "Eif5"
## [41] "Dnajb1" "Srm" "Fam102b" "Eif2s2"
## [45] "P2ry12" "Cebpb" "Calm2" "Set"
## [49] "Ogfrl1" "C5ar1" "Pten" "Phgdh"
## [53] "Bin1" "Eif4a1" "Cct5" "Ppp1r14b"
## [57] "Rab3il1" "Ahsa2" "Rps6ka1" "Jam3"
## [61] "Denr" "St13" "Ifi27" "Serbp1"
## [65] "Irf2bp2" "Tcf7l2" "Mgat2" "Snx18"
## [69] "Npl" "Top1" "Nars" "Fchsd2"
## [73] "Entpd1" "Ptma" "Chsy1" "Ckb"
## [77] "Mrfap1" "Arpc2" "Picalm" "Cct4"
## [81] "Fscn1" "Bag3" "Mcrip1" "Fxr1"
## [85] "Selenos" "Orai3" "Txndc11" "Ywhae"
## [89] "Cdk6" "Ltc4s" "Spsb1" "Pcbp1"
## [93] "Creld2" "Mertk" "Gbp7" "Ranbp1"
## [97] "Nifk" "Scoc" "Klf3" "Magoh"
## [101] "Ran" "Hmgb1" "Npm1" "Tspan7"
## [105] "Bank1" "Tubb5" "Nrip1" "Arhgdia"
## [109] "Plxdc2" "Tmem65" "Spcs2" "Cct2"
## [113] "Pnp" "Ptgs1" "Tmem204" "Pafah1b1"
## [117] "Lsm2" "Ryk" "Cib1" "Gas6"
## [121] "Slc39a10" "Washc4" "Bin2" "Fli1"
## [125] "Sel1l" "Tcea1" "Adgrg1" "Csk"
## [129] "Idh2" "Crebbp" "Ipo7" "Tuba1b"
## [133] "Usp10" "Morf4l2" "Atp2a2" "Clta"
## [137] "Spint1" "Mrps6" "Gcnt1" "G3bp1"
## [141] "Il1a" "Sdf2l1" "Mrps12" "Rab14"
## [145] "Nfam1" "Dynll1" "Vapa" "Cd300a"
## [149] "4833439L19Rik" "Ppa1" "Eif1b" "Ifrd1"
## [153] "B3galt5" "Pa2g4" "Rps19bp1" "Rapgef5"
## [157] "Klhl9" "Eif3c" "Ywhah" "Erh"
## [161] "Pip4k2a" "Pak2" "Sox4" "Ebna1bp2"
## [165] "Cnbp" "Tjp1" "Tmem173" "Alox5ap"
## [169] "Ppp2ca" "Tagln2" "Rbpj" "Mef2a"
## [173] "Rpl36al" "Mak16" "Rangap1" "Eif3j1"
## [177] "Mtus1" "Fam167b" "Rbm8a" "Slc25a5"
## [181] "Eif4e" "Ncl" "Scaf11" "Nudc"
## [185] "Ahsa1" "Hnrnpu" "Smndc1" "Snx6"
## [189] "Fez2" "Mex3c" "Sall1" "Rrp1b"
## [193] "Zfr" "Eif4g2" "Cct8" "Yy1"
## [197] "Gtf2i" "Tmed7" "Mylip" "Cltc"
## [201] "Rexo2" "Comt" "Gmps" "Dctpp1"
## [205] "Gngt2" "C1qbp" "Ubxn4" "Nolc1"
## [209] "Khk"lapply(DEG.list, length)## $LPS.TDT
## [1] 130
##
## $LPS.CTR
## [1] 209GEX.seur <- add_geneset_score(GEX.seur, DEG.list$LPS.TDT, "LPS.TDT.up")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DEG.list$LPS.CTR, "LPS.TDT.dn")## Summarizing datappnew.2.v1 <- VlnPlot(GEX.seur, features = c("LPS.TDT.up","LPS.TDT.dn"),
group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1
ppnew.2.v1a <- VlnPlot(GEX.seur, features = c("LPS.TDT.up","LPS.TDT.dn")[1],
same.y.lims = T,
group.by = "cnt", cols = color.cnt,
ncol = 1, pt.size = 0, raster = F) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P30.PBS.CTR","P30.PBS.MIG"),
c("P30.PBS.MIG","P30.PBS.TDT"),
c("P30.PBS.CTR","P30.PBS.TDT"),
c("P30.LPS.CTR","P30.LPS.MIG"),
c("P30.LPS.MIG","P30.LPS.TDT"),
c("P30.LPS.CTR","P30.LPS.TDT"),
c("P30.PBS.TDT","P30.LPS.TDT")),
label.y = c(0.55,0.75,0.95,0.55,0.75,0.95,1.1),
size=3
)
ppnew.2.v1a
ppnew.2.v1b <- VlnPlot(GEX.seur, features = c("LPS.TDT.up","LPS.TDT.dn")[2],
same.y.lims = T,
group.by = "cnt", cols = color.cnt,
ncol = 1, pt.size = 0, raster = F) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & ylim(c(-0.23,0.57)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P30.PBS.CTR","P30.PBS.MIG"),
c("P30.PBS.MIG","P30.PBS.TDT"),
c("P30.PBS.CTR","P30.PBS.TDT"),
c("P30.LPS.CTR","P30.LPS.MIG"),
c("P30.LPS.MIG","P30.LPS.TDT"),
c("P30.LPS.CTR","P30.LPS.TDT")#,
#c("P30.PBS.TDT","P30.LPS.TDT")
),
label.y = c(0.25,0.3,0.4,0.35,0.4,0.5,0.55),
size=3
)
ppnew.2.v1b
DotPlot(GEX.seur, features = c("LPS.TDT.dn","LPS.TDT.up"), group.by = "cnt", cols = c("midnightblue","darkorange1"), scale = F, scale.by = "size", scale.min = 1, scale.max = 1, dot.scale = 12) +
coord_flip() + labs(title="mean signature score") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) 
GEX.seur@meta.data[,c("cnt","LPS.TDT.dn","LPS.TDT.up")] %>%
reshape2::melt("cnt") %>%
ggplot(aes(x = cnt, y = variable, color=value)) +
geom_point(size=12) +
scale_color_gradientn(colors = c("midnightblue","darkorange1")) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
new clusters
GEX.seur$new_clusters <- as.character(GEX.seur$seurat_clusters)
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(1)] <- "C1"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(2,8,7,14)] <- "C2"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(9)] <- "C3"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(12)] <- "C4"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(3)] <- "C5"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(5)] <- "C6"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(0,4)] <- "C7"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(6)] <- "C8"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(11)] <- "C9"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(13)] <- "C10"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(10)] <- "C11"
GEX.seur$new_clusters <- factor(GEX.seur$new_clusters,
levels = paste0("C",1:11))head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGGTCGACA-1 Spp1tdt 5001 2062 2.699460 12.437512
## AAACCCAAGGTTCTTG-1 Spp1tdt 6501 2457 2.553453 12.305799
## AAACCCACAAACTCGT-1 Spp1tdt 4188 1841 2.960840 8.572111
## AAACCCACACGCGGTT-1 Spp1tdt 4524 1913 1.215738 9.858532
## AAACCCACATCCGTGG-1 Spp1tdt 2337 1225 3.979461 7.659392
## AAACCCAGTATCGAAA-1 Spp1tdt 4318 1955 2.292728 11.741547
## FB.info S.Score G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGGTCGACA-1 P30.LPS.CTR 0.0061184939 -0.06523366 S 0
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT1 0.0002076412 -0.07102627 S 0
## AAACCCACAAACTCGT-1 P30.PBS.MIG 0.0192137320 0.06031440 G2M 7
## AAACCCACACGCGGTT-1 P30.LPS.TDT2 -0.0535299003 -0.05464656 G1 0
## AAACCCACATCCGTGG-1 P30.PBS.TDT 0.0489756368 -0.04211671 S 2
## AAACCCAGTATCGAAA-1 P30.LPS.MIG -0.0314645626 0.02183367 G2M 0
## seurat_clusters pANN_0.25_0.005_1191 DoubletFinder0.05
## AAACCCAAGGTCGACA-1 0 0.000000000 Singlet
## AAACCCAAGGTTCTTG-1 0 0.044025157 Singlet
## AAACCCACAAACTCGT-1 7 0.006289308 Singlet
## AAACCCACACGCGGTT-1 0 0.012578616 Singlet
## AAACCCACATCCGTGG-1 2 0.000000000 Singlet
## AAACCCAGTATCGAAA-1 0 0.006289308 Singlet
## pANN_0.25_0.005_2383 DoubletFinder0.1 sort_clusters age
## AAACCCAAGGTCGACA-1 0.000000000 Singlet 0 P30.LPS
## AAACCCAAGGTTCTTG-1 0.044025157 Singlet 0 P30.LPS
## AAACCCACAAACTCGT-1 0.006289308 Singlet 7 P30.PBS
## AAACCCACACGCGGTT-1 0.012578616 Singlet 0 P30.LPS
## AAACCCACATCCGTGG-1 0.000000000 Singlet 2 P30.PBS
## AAACCCAGTATCGAAA-1 0.006289308 Singlet 0 P30.LPS
## cnt preAnno RNA_snn_res.1 RNA_snn_res.0.8
## AAACCCAAGGTCGACA-1 P30.LPS.CTR Microglia 2 2
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT Microglia 3 4
## AAACCCACAAACTCGT-1 P30.PBS.MIG Microglia 7 0
## AAACCCACACGCGGTT-1 P30.LPS.TDT Microglia 2 2
## AAACCCACATCCGTGG-1 P30.PBS.TDT Microglia 0 0
## AAACCCAGTATCGAAA-1 P30.LPS.MIG Microglia 2 2
## RNA_snn_res.0.6 FB.new RNA_snn_res.1.2 RNA_snn_res.1.3
## AAACCCAAGGTCGACA-1 0 P30.LPS.CTR 2 2
## AAACCCAAGGTTCTTG-1 0 P30.LPS.TDT1 3 4
## AAACCCACAAACTCGT-1 1 P30.PBS.MIG 7 7
## AAACCCACACGCGGTT-1 0 P30.LPS.TDT2 3 4
## AAACCCACATCCGTGG-1 1 P30.PBS.TDT 1 1
## AAACCCAGTATCGAAA-1 0 P30.LPS.MIG 2 2
## RNA_snn_res.1.4 LPS.TDT.up LPS.TDT.dn new_clusters age1
## AAACCCAAGGTCGACA-1 2 0.2671683 0.24670044 C7 P30.LPS
## AAACCCAAGGTTCTTG-1 2 0.2190005 0.21501295 C7 P30.LPS
## AAACCCACAAACTCGT-1 7 0.3136800 0.02545367 C2 P30.PBS
## AAACCCACACGCGGTT-1 2 0.2974201 0.09338095 C7 P30.LPS
## AAACCCACATCCGTGG-1 1 0.2758054 -0.03406151 C2 P30.PBS
## AAACCCAGTATCGAAA-1 2 0.2300058 0.11176488 C7 P30.LPS
## cluster_cnt DAM.sig_top50 DAM.sig_top100 DAM.sig_top250
## AAACCCAAGGTCGACA-1 C7_P30.LPS.CTR 0.12647123 0.09610184 0.09546331
## AAACCCAAGGTTCTTG-1 C7_P30.LPS.TDT 0.04136910 0.02314892 0.06780191
## AAACCCACAAACTCGT-1 C2_P30.PBS.MIG -0.03826143 -0.02304032 -0.01821938
## AAACCCACACGCGGTT-1 C7_P30.LPS.TDT 0.17479143 0.26609309 0.19238189
## AAACCCACATCCGTGG-1 C2_P30.PBS.TDT 0.06229684 -0.06607385 -0.05013960
## AAACCCAGTATCGAAA-1 C7_P30.LPS.MIG 0.19221428 0.07223828 0.07965503
## DAM.sig_top500 LPS.sig_score
## AAACCCAAGGTCGACA-1 0.2497200 0.5988063
## AAACCCAAGGTTCTTG-1 0.2634468 0.6298800
## AAACCCACAAACTCGT-1 0.1498645 0.2308747
## AAACCCACACGCGGTT-1 0.3066293 0.6339013
## AAACCCACATCCGTGG-1 0.1023302 0.1488521
## AAACCCAGTATCGAAA-1 0.2186862 0.4894010#
GEX.seur$age1 <- factor(as.character(GEX.seur$age),
levels = c("P30.PBS","P30.LPS"))
#
GEX.seur$FB.new <- factor(as.character(GEX.seur$FB.info),
levels = rev(levels(GEX.seur$FB.info))[c(8:10,4,5,7,6)])
color.FBnew <- color.FB2[c(3:1,7,6,4,5)]head(GEX.seur$cnt)## AAACCCAAGGTCGACA-1 AAACCCAAGGTTCTTG-1 AAACCCACAAACTCGT-1 AAACCCACACGCGGTT-1
## P30.LPS.CTR P30.LPS.TDT P30.PBS.MIG P30.LPS.TDT
## AAACCCACATCCGTGG-1 AAACCCAGTATCGAAA-1
## P30.PBS.TDT P30.LPS.MIG
## 6 Levels: P30.PBS.CTR P30.PBS.MIG P30.PBS.TDT P30.LPS.CTR ... P30.LPS.TDTDotPlot(GEX.seur, features = c(markers.mig,"Ifit1","Ifit2"), group.by = "new_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
new markers
# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "new_clusters"
#GEX.markers.new <- FindAllMarkers(GEX.seur, only.pos = TRUE,
# min.pct = 0.05,
# test.use = "MAST",
# logfc.threshold = 0.25)
GEX.markers.new <- read.table("sc10x_LYN.0921.marker.subset2_new.csv", header = TRUE, sep = ",")
GEX.markers.new %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)## # A tibble: 88 x 7
## # Groups: cluster [11]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 0 1.05 0.999 0.973 0 C1 P2ry12
## 2 0 0.965 0.998 0.921 0 C1 Gpr34
## 3 0 0.812 0.979 0.897 0 C1 Tgfbr1
## 4 2.09e-269 0.751 0.966 0.857 3.73e-265 C1 Arhgap5
## 5 1.65e-265 0.731 0.973 0.903 2.96e-261 C1 Maf
## 6 4.40e-250 0.748 0.925 0.764 7.88e-246 C1 Rhob
## 7 9.96e-142 0.717 0.511 0.276 1.78e-137 C1 Slc40a1
## 8 3.05e-121 0.717 0.566 0.329 5.46e-117 C1 Sox4
## 9 0 1.03 0.993 0.91 0 C2 Gpr34
## 10 0 1.01 1 0.968 0 C2 P2ry12
## # ... with 78 more rowsmarkers.new <- GEX.markers.new
markers.new$cluster <- factor(as.character(GEX.markers.new$cluster),
levels = levels(GEX.seur$new_clusters))
markers.new_t60 <- (markers.new %>% group_by(cluster) %>%
filter(pct.1>0.1 & gene %in% grep("Rps|Rpl|Gm|Rik$|Xist|Tsix",markers.new$gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
ungroup() %>%
#arrange(desc(avg_log2FC*pct.1),gene) %>%
arrange(desc(pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$geneDotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[1:64])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[65:128])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[129:192])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[193:256])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[257:320])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[321:384])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[385:436])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
plot1
umap in clusters
pp.umap.1a <- DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters", raster = F, pt.size = 0.3) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.new", cols = color.FBnew, raster = F, pt.size = 0.3)
pp.umap.1a
pp.umap.1b <- DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "new_clusters", repel = F, raster = F, pt.size = 0.3) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt, raster = F, pt.size = 0.3)
pp.umap.1b
plot2
pp.umap.2a <- DimPlot(GEX.seur, label = F, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", raster = F, pt.size = 0.3,
ncol = 3, cols = color.cnt)
pp.umap.2a
pp.umap.2b <- DimPlot(GEX.seur, label = F,repel = F, reduction = 'umap', group.by = "cnt", split.by = "age",
raster = F, pt.size = 0.3,
ncol = 3, cols = color.cnt)
pp.umap.2b
pp.umap.2c <- DimPlot(GEX.seur, label = F, repel = F, reduction = 'umap', group.by = "new_clusters", split.by = "cnt", raster = F, pt.size = 0.3,
ncol = 3)
pp.umap.2c
pp.umap.2d <- DimPlot(GEX.seur, label = F,repel = F, reduction = 'umap', group.by = "new_clusters", split.by = "age",
raster = F, pt.size = 0.3,
ncol = 3)
pp.umap.2d
plot3
composition
pp.comp.3a <- cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$cnt,
clusters=GEX.seur$new_clusters),
main = "Cell Count",
gaps_row = c(3),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$cnt,
clusters=GEX.seur$new_clusters)),
main = "Cell Ratio",
gaps_row = c(3),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
pp.comp.3a
stat_sort <- GEX.seur@meta.data[,c("cnt","new_clusters")]stat_sort.s <- stat_sort %>%
group_by(cnt,new_clusters) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup()## `summarise()` has grouped output by 'cnt'. You can override using the `.groups`
## argument.dim(stat_sort.s)## [1] 55 4#stat_sort.spp.comp.3b <- stat_sort.s %>%
ggplot(aes(x = new_clusters, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.FBnew, name="") +
labs(title="Cell Composition", y="Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x = new_clusters, y = 100*prop, color=cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=1.75,
stroke=0.15, show.legend = F)
pp.comp.3b
plot4
features4 <- c("Saa3","Il1a","Ccl3","Ccl4",
"Ccl5","Il1b","Tnf","Ccl2",
"Cxcl10")pp.feat.4 <- FeaturePlot(GEX.seur,
features = features4, raster = T, pt.size = 3.5,
ncol = 4)
pp.feat.4
plot5
final.markers <- c("Apoe","Ms4a6d","Ms4a6c","Ctss",
"Ctsc","Ctsd","Ch25h","Ccl2",
"Ccl3","Ccl4","Ccl5","Cxcl10",
"Il1a","Il1b","Tnf","Saa3")
final.markers## [1] "Apoe" "Ms4a6d" "Ms4a6c" "Ctss" "Ctsc" "Ctsd" "Ch25h" "Ccl2"
## [9] "Ccl3" "Ccl4" "Ccl5" "Cxcl10" "Il1a" "Il1b" "Tnf" "Saa3"length(final.markers)## [1] 16sum(final.markers %in% markers.new$gene)## [1] 16final.markers.sort <- (markers.new %>% group_by(cluster) %>%
filter(gene %in% final.markers) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster, p_val_adj))$gene
final.markers.sort## [1] "Ctsd" "Il1a" "Saa3" "Ms4a6c" "Ch25h" "Ms4a6d" "Ctsc" "Ccl5"
## [9] "Ccl3" "Tnf" "Il1b" "Cxcl10" "Ccl2" "Ccl4" "Apoe" "Ctss"final.markers.new <- final.markers.sort
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)pp.dot.5a1 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters")
pp.dot.5a1
pp.dot.5a2 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters") +
scale_color_gradientn(colours = rev(mapal))
pp.dot.5a2
pp.dot.5b1 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "cnt"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Conditions")
pp.dot.5b1
pp.dot.5b2 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "cnt"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Conditions") +
scale_color_gradientn(colours = rev(mapal))
pp.dot.5b2
levels.cluster_cnt <- as.vector(unlist(sapply(levels(GEX.seur$new_clusters),function(x){
paste0(x,"_",levels(GEX.seur$cnt))
})))
levels.cluster_cnt## [1] "C1_P30.PBS.CTR" "C1_P30.PBS.MIG" "C1_P30.PBS.TDT" "C1_P30.LPS.CTR"
## [5] "C1_P30.LPS.MIG" "C1_P30.LPS.TDT" "C2_P30.PBS.CTR" "C2_P30.PBS.MIG"
## [9] "C2_P30.PBS.TDT" "C2_P30.LPS.CTR" "C2_P30.LPS.MIG" "C2_P30.LPS.TDT"
## [13] "C3_P30.PBS.CTR" "C3_P30.PBS.MIG" "C3_P30.PBS.TDT" "C3_P30.LPS.CTR"
## [17] "C3_P30.LPS.MIG" "C3_P30.LPS.TDT" "C4_P30.PBS.CTR" "C4_P30.PBS.MIG"
## [21] "C4_P30.PBS.TDT" "C4_P30.LPS.CTR" "C4_P30.LPS.MIG" "C4_P30.LPS.TDT"
## [25] "C5_P30.PBS.CTR" "C5_P30.PBS.MIG" "C5_P30.PBS.TDT" "C5_P30.LPS.CTR"
## [29] "C5_P30.LPS.MIG" "C5_P30.LPS.TDT" "C6_P30.PBS.CTR" "C6_P30.PBS.MIG"
## [33] "C6_P30.PBS.TDT" "C6_P30.LPS.CTR" "C6_P30.LPS.MIG" "C6_P30.LPS.TDT"
## [37] "C7_P30.PBS.CTR" "C7_P30.PBS.MIG" "C7_P30.PBS.TDT" "C7_P30.LPS.CTR"
## [41] "C7_P30.LPS.MIG" "C7_P30.LPS.TDT" "C8_P30.PBS.CTR" "C8_P30.PBS.MIG"
## [45] "C8_P30.PBS.TDT" "C8_P30.LPS.CTR" "C8_P30.LPS.MIG" "C8_P30.LPS.TDT"
## [49] "C9_P30.PBS.CTR" "C9_P30.PBS.MIG" "C9_P30.PBS.TDT" "C9_P30.LPS.CTR"
## [53] "C9_P30.LPS.MIG" "C9_P30.LPS.TDT" "C10_P30.PBS.CTR" "C10_P30.PBS.MIG"
## [57] "C10_P30.PBS.TDT" "C10_P30.LPS.CTR" "C10_P30.LPS.MIG" "C10_P30.LPS.TDT"
## [61] "C11_P30.PBS.CTR" "C11_P30.PBS.MIG" "C11_P30.PBS.TDT" "C11_P30.LPS.CTR"
## [65] "C11_P30.LPS.MIG" "C11_P30.LPS.TDT"GEX.seur$cluster_cnt <- factor(paste0(as.character(GEX.seur$new_clusters),
"_",
as.character(GEX.seur$cnt)),
levels = levels.cluster_cnt)pp.dot.5c1a <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)),
features = final.markers.new, group.by = "cluster_cnt"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions")
pp.dot.5c1a
pp.dot.5c2a <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)),
features = final.markers.new, group.by = "cluster_cnt"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions") +
scale_color_gradientn(colours = rev(mapal))
pp.dot.5c2a
pp.dot.5c1b <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)),
features = final.markers.new, group.by = "cluster_cnt"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions")
pp.dot.5c1b
pp.dot.5c2b <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)),
features = final.markers.new, group.by = "cluster_cnt"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions") +
scale_color_gradientn(colours = rev(mapal))
pp.dot.5c2b
plot6
violin using same markers in plot5
pp.violin.6a1 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters",
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6a1
pp.violin.6a2 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters",
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6a2
pp.violin.6b1 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6b1
pp.violin.6b2 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6b2
pp.violin.6b3 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0,7.2)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
c("P30.LPS.MIG","P30.LPS.TDT"),
c("P30.LPS.CTR","P30.LPS.TDT"),
c("P30.PBS.CTR","P30.PBS.MIG"),
c("P30.PBS.MIG","P30.PBS.TDT"),
c("P30.PBS.CTR","P30.PBS.TDT"),
c("P30.LPS.TDT","P30.PBS.TDT")),
label.y = c(6.3,6.8,7.3,6.3,6.8,7.3,7.8)-1,
size=1.75
)
pp.violin.6b3
pp.violin.6c1a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)),
features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.FBnew[1:3],4),
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6c1a
pp.violin.6c1b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)),
features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6c1b
pp.violin.6c2a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)),
features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.FBnew[1:3],4),
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6c2a
pp.violin.6c2b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)),
features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6c2b
contrast.a <- list(c("P30.PBS.CTR","P30.PBS.MIG"),
c("P30.PBS.MIG","P30.PBS.TDT"),
c("P30.PBS.CTR","P30.PBS.TDT"))
plot.contrast.a <- list()
for(XXX in paste0("C",1:3)){
plot.contrast.a <- c(plot.contrast.a,
lapply(contrast.a, function(x){
paste0(XXX,
"_",
x)
}))
}
plot.labely.a <- rep(c(6.3,6.8,7.3),3)
#
length(plot.contrast.a)## [1] 9plot.contrast.a[1:6]## [[1]]
## [1] "C1_P30.PBS.CTR" "C1_P30.PBS.MIG"
##
## [[2]]
## [1] "C1_P30.PBS.MIG" "C1_P30.PBS.TDT"
##
## [[3]]
## [1] "C1_P30.PBS.CTR" "C1_P30.PBS.TDT"
##
## [[4]]
## [1] "C2_P30.PBS.CTR" "C2_P30.PBS.MIG"
##
## [[5]]
## [1] "C2_P30.PBS.MIG" "C2_P30.PBS.TDT"
##
## [[6]]
## [1] "C2_P30.PBS.CTR" "C2_P30.PBS.TDT"contrast.b <- list(c("P30.LPS.CTR","P30.LPS.MIG"),
c("P30.LPS.MIG","P30.LPS.TDT"),
c("P30.LPS.CTR","P30.LPS.TDT"))
plot.contrast.b <- list()
for(XXX in paste0("C",4:11)){
plot.contrast.b <- c(plot.contrast.b,
lapply(contrast.b, function(x){
paste0(XXX,
"_",
x)
}))
}
plot.labely.b <- rep(c(6.3,6.8,7.3),8)
#
length(plot.contrast.b)## [1] 24plot.contrast.b[1:6]## [[1]]
## [1] "C4_P30.LPS.CTR" "C4_P30.LPS.MIG"
##
## [[2]]
## [1] "C4_P30.LPS.MIG" "C4_P30.LPS.TDT"
##
## [[3]]
## [1] "C4_P30.LPS.CTR" "C4_P30.LPS.TDT"
##
## [[4]]
## [1] "C5_P30.LPS.CTR" "C5_P30.LPS.MIG"
##
## [[5]]
## [1] "C5_P30.LPS.MIG" "C5_P30.LPS.TDT"
##
## [[6]]
## [1] "C5_P30.LPS.CTR" "C5_P30.LPS.TDT"pp.violin.6c3a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)),
features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.FBnew[1:3],4),
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)& ylim(c(0,6.8)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = plot.contrast.a,
label.y = plot.labely.a-1,
size=1.5
)
pp.violin.6c3a
pp.violin.6c3b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)),
features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)& ylim(c(0,6.8)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = plot.contrast.b,
label.y = plot.labely.b-1,
size=1.75
)
pp.violin.6c3b
plot7
count DEG number, using different cutoffs
cut.padj = 0.05
cut.log2FC = 0.1
cut.pct1 = 0.05
stat.a1 <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat.a1## contrast cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR 209
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT 130
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG 187
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT 164
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR 28
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT 35
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG 21
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT 5#df.DEGs.app.stat.DEG[["a1"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt[c(4:6,1:3)], name="") +
labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["a1"]]
cut.padj = 0.05
cut.log2FC = 0.2
cut.pct1 = 0.05
stat.a2 <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat.a2## contrast cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR 30
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT 53
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG 44
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT 94
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR 9
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT 17
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG 9
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT 4pp.stat.DEG[["a2"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt[c(4:6,1:3)], name="") +
labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["a2"]]
cut.padj = 0.05
cut.log2FC = 0.3
cut.pct1 = 0.05
stat.a3 <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat.a3## contrast cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR 8
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT 29
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG 17
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT 51
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR 2
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT 7
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG 3
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT 1pp.stat.DEG[["a3"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt[c(4:6,1:3)], name="") +
labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["a3"]]
DEG list
cut.padj = 0.05
cut.log2FC = 0.1
cut.pct1 = 0.05
stat.a1.df <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
#summarise(up.DEGs = n()) %>%
as.data.frame()
dim(stat.a1.df)## [1] 779 8stat.a1.df[1:6,]## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 7.291059e-54 0.5042723 0.174 0.033 1.305100e-49 P30.PBS.CTR tdTomato-head
## 2 5.510551e-20 0.1877754 0.998 0.990 9.863886e-16 P30.PBS.CTR Fcrls
## 3 7.112810e-19 0.3197571 0.447 0.310 1.273193e-14 P30.PBS.CTR Lrba
## 4 2.383831e-13 0.1595905 0.999 0.999 4.267058e-09 P30.PBS.CTR P2ry12
## 5 7.947888e-13 0.1056330 1.000 1.000 1.422672e-08 P30.PBS.CTR Actb
## 6 1.076088e-12 0.2209830 0.892 0.837 1.926198e-08 P30.PBS.CTR Cd164
## contrast
## 1 P30.PBS.TDTvsCTR
## 2 P30.PBS.TDTvsCTR
## 3 P30.PBS.TDTvsCTR
## 4 P30.PBS.TDTvsCTR
## 5 P30.PBS.TDTvsCTR
## 6 P30.PBS.TDTvsCTRcut.padj = 0.05
cut.log2FC = 0.2
cut.pct1 = 0.05
stat.a2.df <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
#summarise(up.DEGs = n()) %>%
as.data.frame()
dim(stat.a2.df)## [1] 260 8stat.a2.df[1:6,]## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 7.291059e-54 0.5042723 0.174 0.033 1.305100e-49 P30.PBS.CTR tdTomato-head
## 2 7.112810e-19 0.3197571 0.447 0.310 1.273193e-14 P30.PBS.CTR Lrba
## 3 1.076088e-12 0.2209830 0.892 0.837 1.926198e-08 P30.PBS.CTR Cd164
## 4 2.513301e-12 0.2151293 0.892 0.828 4.498809e-08 P30.PBS.CTR Pmepa1
## 5 3.785970e-12 0.2683232 0.287 0.192 6.776886e-08 P30.PBS.CTR Bank1
## 6 2.500566e-11 0.2553055 0.137 0.071 4.476013e-07 P30.PBS.CTR Upk1b
## contrast
## 1 P30.PBS.TDTvsCTR
## 2 P30.PBS.TDTvsCTR
## 3 P30.PBS.TDTvsCTR
## 4 P30.PBS.TDTvsCTR
## 5 P30.PBS.TDTvsCTR
## 6 P30.PBS.TDTvsCTRcut.padj = 0.05
cut.log2FC = 0.3
cut.pct1 = 0.05
stat.a3.df <- df.DEGs.a %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(contrast,cluster) %>%
#summarise(up.DEGs = n()) %>%
as.data.frame()
dim(stat.a3.df)## [1] 118 8stat.a3.df[1:6,]## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 7.291059e-54 0.5042723 0.174 0.033 1.305100e-49 P30.PBS.CTR tdTomato-head
## 2 7.112810e-19 0.3197571 0.447 0.310 1.273193e-14 P30.PBS.CTR Lrba
## 3 3.177818e-55 0.6194995 0.827 0.722 5.688295e-51 P30.PBS.TDT H2-D1
## 4 1.418235e-47 0.5638234 0.819 0.743 2.538641e-43 P30.PBS.TDT H2-K1
## 5 1.868002e-28 0.7477478 0.966 0.957 3.343723e-24 P30.PBS.TDT Apoe
## 6 1.724448e-20 0.5015367 0.495 0.400 3.086762e-16 P30.PBS.TDT Lyz2
## contrast
## 1 P30.PBS.TDTvsCTR
## 2 P30.PBS.TDTvsCTR
## 3 P30.PBS.TDTvsCTR
## 4 P30.PBS.TDTvsCTR
## 5 P30.PBS.TDTvsCTR
## 6 P30.PBS.TDTvsCTR#
stat.a1.list <- list(P30.LPS.TDTvsCTR.CTRup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.CTR"))$gene,
P30.LPS.TDTvsCTR.TDTup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.TDT"))$gene,
P30.LPS.TDTvsMIG.MIGup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.MIG"))$gene,
P30.LPS.TDTvsMIG.TDTup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.TDT"))$gene,
P30.PBS.TDTvsCTR.CTRup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.CTR"))$gene,
P30.PBS.TDTvsCTR.TDTup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.TDT"))$gene,
P30.PBS.TDTvsMIG.MIGup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.MIG"))$gene,
P30.PBS.TDTvsMIG.TDTup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.TDT"))$gene)
stat.a2.list <- list(P30.LPS.TDTvsCTR.CTRup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.CTR"))$gene,
P30.LPS.TDTvsCTR.TDTup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.TDT"))$gene,
P30.LPS.TDTvsMIG.MIGup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.MIG"))$gene,
P30.LPS.TDTvsMIG.TDTup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.TDT"))$gene,
P30.PBS.TDTvsCTR.CTRup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.CTR"))$gene,
P30.PBS.TDTvsCTR.TDTup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.TDT"))$gene,
P30.PBS.TDTvsMIG.MIGup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.MIG"))$gene,
P30.PBS.TDTvsMIG.TDTup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.TDT"))$gene)
stat.a3.list <- list(P30.LPS.TDTvsCTR.CTRup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.CTR"))$gene,
P30.LPS.TDTvsCTR.TDTup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.TDT"))$gene,
P30.LPS.TDTvsMIG.MIGup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.MIG"))$gene,
P30.LPS.TDTvsMIG.TDTup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.TDT"))$gene,
P30.PBS.TDTvsCTR.CTRup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.CTR"))$gene,
P30.PBS.TDTvsCTR.TDTup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.TDT"))$gene,
P30.PBS.TDTvsMIG.MIGup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.MIG"))$gene,
P30.PBS.TDTvsMIG.TDTup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.TDT"))$gene)lapply(stat.a1.list,length)## $P30.LPS.TDTvsCTR.CTRup
## [1] 209
##
## $P30.LPS.TDTvsCTR.TDTup
## [1] 130
##
## $P30.LPS.TDTvsMIG.MIGup
## [1] 187
##
## $P30.LPS.TDTvsMIG.TDTup
## [1] 164
##
## $P30.PBS.TDTvsCTR.CTRup
## [1] 28
##
## $P30.PBS.TDTvsCTR.TDTup
## [1] 35
##
## $P30.PBS.TDTvsMIG.MIGup
## [1] 21
##
## $P30.PBS.TDTvsMIG.TDTup
## [1] 5lapply(stat.a2.list,length)## $P30.LPS.TDTvsCTR.CTRup
## [1] 30
##
## $P30.LPS.TDTvsCTR.TDTup
## [1] 53
##
## $P30.LPS.TDTvsMIG.MIGup
## [1] 44
##
## $P30.LPS.TDTvsMIG.TDTup
## [1] 94
##
## $P30.PBS.TDTvsCTR.CTRup
## [1] 9
##
## $P30.PBS.TDTvsCTR.TDTup
## [1] 17
##
## $P30.PBS.TDTvsMIG.MIGup
## [1] 9
##
## $P30.PBS.TDTvsMIG.TDTup
## [1] 4lapply(stat.a3.list,length)## $P30.LPS.TDTvsCTR.CTRup
## [1] 8
##
## $P30.LPS.TDTvsCTR.TDTup
## [1] 29
##
## $P30.LPS.TDTvsMIG.MIGup
## [1] 17
##
## $P30.LPS.TDTvsMIG.TDTup
## [1] 51
##
## $P30.PBS.TDTvsCTR.CTRup
## [1] 2
##
## $P30.PBS.TDTvsCTR.TDTup
## [1] 7
##
## $P30.PBS.TDTvsMIG.MIGup
## [1] 3
##
## $P30.PBS.TDTvsMIG.TDTup
## [1] 1load subset1 DEGs.a list
#
stat.s1.list <- list(P06.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsCTR.CTRup.txt"))),
P06.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsCTR.TDTup.txt"))),
P06.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsMIG.MIGup.txt"))),
P06.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsMIG.TDTup.txt"))),
P30.PBS.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsCTR.CTRup.txt"))),
P30.PBS.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsCTR.TDTup.txt"))),
P30.PBS.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsMIG.MIGup.txt"))),
P30.PBS.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsMIG.TDTup.txt"))))
#
stat.s2.list <- list(P06.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsCTR.CTRup.txt"))),
P06.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsCTR.TDTup.txt"))),
P06.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsMIG.MIGup.txt"))),
P06.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsMIG.TDTup.txt"))),
P30.PBS.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsCTR.CTRup.txt"))),
P30.PBS.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsCTR.TDTup.txt"))),
P30.PBS.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsMIG.MIGup.txt"))),
P30.PBS.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsMIG.TDTup.txt"))))
#
stat.s3.list <- list(P06.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsCTR.CTRup.txt"))),
P06.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsCTR.TDTup.txt"))),
P06.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsMIG.MIGup.txt"))),
P06.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsMIG.TDTup.txt"))),
P30.PBS.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsCTR.CTRup.txt"))),
P30.PBS.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsCTR.TDTup.txt"))),
P30.PBS.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsMIG.MIGup.txt"))),
P30.PBS.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsMIG.TDTup.txt"))))lapply(stat.s1.list,length)## $P06.TDTvsCTR.CTRup
## [1] 174
##
## $P06.TDTvsCTR.TDTup
## [1] 84
##
## $P06.TDTvsMIG.MIGup
## [1] 70
##
## $P06.TDTvsMIG.TDTup
## [1] 61
##
## $P30.PBS.TDTvsCTR.CTRup
## [1] 28
##
## $P30.PBS.TDTvsCTR.TDTup
## [1] 35
##
## $P30.PBS.TDTvsMIG.MIGup
## [1] 21
##
## $P30.PBS.TDTvsMIG.TDTup
## [1] 5lapply(stat.s2.list,length)## $P06.TDTvsCTR.CTRup
## [1] 62
##
## $P06.TDTvsCTR.TDTup
## [1] 52
##
## $P06.TDTvsMIG.MIGup
## [1] 11
##
## $P06.TDTvsMIG.TDTup
## [1] 31
##
## $P30.PBS.TDTvsCTR.CTRup
## [1] 8
##
## $P30.PBS.TDTvsCTR.TDTup
## [1] 17
##
## $P30.PBS.TDTvsMIG.MIGup
## [1] 9
##
## $P30.PBS.TDTvsMIG.TDTup
## [1] 4lapply(stat.s3.list,length)## $P06.TDTvsCTR.CTRup
## [1] 14
##
## $P06.TDTvsCTR.TDTup
## [1] 18
##
## $P06.TDTvsMIG.MIGup
## [1] 3
##
## $P06.TDTvsMIG.TDTup
## [1] 14
##
## $P30.PBS.TDTvsCTR.CTRup
## [1] 2
##
## $P30.PBS.TDTvsCTR.TDTup
## [1] 7
##
## $P30.PBS.TDTvsMIG.MIGup
## [1] 3
##
## $P30.PBS.TDTvsMIG.TDTup
## [1] 1plot8
DEG volcano, set same x/y-axis
pp.DEGs.a$P30.LPS.TDTvsCTR <- DEGs.a.combine$P30.LPS.TDTvsCTR %>%
mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTR"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.LPS, TDT vs CTR") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.LPS.TDTvsCTR
pp.DEGs.a$P30.LPS.TDTvsMIG <- DEGs.a.combine$P30.LPS.TDTvsMIG %>%
mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("MIG"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.LPS, TDT vs MIG") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.LPS.TDTvsMIG
pp.DEGs.a$P30.PBS.TDTvsCTR <- DEGs.a.combine$P30.PBS.TDTvsCTR %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTR"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.PBS, TDT vs CTR") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.PBS.TDTvsCTR
pp.DEGs.a$P30.PBS.TDTvsMIG <- DEGs.a.combine$P30.PBS.TDTvsMIG %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("MIG"="Skyblue",
"TDT"="pink",
"None"="grey")) +
theme_classic() + labs(title="P30.PBS, TDT vs MIG") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.PBS.TDTvsMIG
plot9
‘LPS_Tdt vs LPS_MIG’, how many DEGs overlap with ‘P06 Tdt vs P06 CTRL’ and ‘P06 Tdt vs MIG’ ?
library(UpSetR)## Warning: package 'UpSetR' was built under R version 4.0.5upset(fromList(c(stat.a1.list[3:4],
stat.s1.list[1:4])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text("pct.a>0.05, padj<0.05, |log2FC|>0.1",x = 0.75, y=0.95, gp=grid::gpar(fontsize=20))
upset(fromList(c(stat.a1.list[3:4],
stat.s1.list[1:2])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text("pct.a>0.05, padj<0.05, |log2FC|>0.1",x = 0.75, y=0.95, gp=grid::gpar(fontsize=20))
upset(fromList(c(stat.a1.list[3:4],
stat.s1.list[3:4])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text("pct.a>0.05, padj<0.05, |log2FC|>0.1",x = 0.75, y=0.95, gp=grid::gpar(fontsize=20))
overlap list
plot10
cut.padj = 0.05
cut.log2FC = 0.3
cut.pct1 = 0.05
df.DEGs.b %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(condition,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()## condition cluster up.DEGs
## 1 P30.CTR P30.PBS 425
## 2 P30.CTR P30.LPS 512
## 3 P30.MIG P30.PBS 321
## 4 P30.MIG P30.LPS 391
## 5 P30.TDT P30.PBS 388
## 6 P30.TDT P30.LPS 423pp.stat.DEG[["b"]] <- df.DEGs.b %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(condition,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=condition, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
pp.stat.DEG[["b"]]
sig.LPS <- df.DEGs.b %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
filter(condition == "P30.CTR" & cluster=="P30.LPS")
dim(sig.LPS)## [1] 512 8sig.LPS[1:6,]## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene condition
## Ccl12 0 3.787020 0.993 0.731 0 P30.LPS Ccl12 P30.CTR
## Ifitm3 0 2.120287 0.766 0.146 0 P30.LPS Ifitm3 P30.CTR
## Slfn2 0 1.709130 0.842 0.190 0 P30.LPS Slfn2 P30.CTR
## Ifi30 0 1.647601 0.966 0.558 0 P30.LPS Ifi30 P30.CTR
## Hsp90ab1 0 1.617709 1.000 0.937 0 P30.LPS Hsp90ab1 P30.CTR
## Srgn 0 1.541016 0.942 0.499 0 P30.LPS Srgn P30.CTRGEX.seur <- add_geneset_score(GEX.seur, sig.LPS$gene, "LPS.sig_score")## Summarizing datapp.feat.10a1 <- FeaturePlot(GEX.seur,features = c("LPS.sig_score"), ncol = 1,
raster = T, pt.size = 3.5#,keep.scale = "all"
)
pp.feat.10a1
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
pp.feat.10a2 <- FeaturePlot(GEX.seur,features = c("LPS.sig_score"), ncol = 1,
raster = T, pt.size = 3.5#,keep.scale = "all"
) & scale_color_gradientn(colors = rev(mapal))
pp.feat.10a2
pp.feat.10b1 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"),
group.by = "cnt", cols = color.cnt,
ncol = 1, pt.size = 0, raster = F) + geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02) + NoLegend()
pp.feat.10b1
pp.feat.10b2 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"),
same.y.lims = T,
group.by = "cnt", cols = color.cnt,
ncol = 1, pt.size = 0, raster = F) + NoLegend() &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & ylim(c(0,1.1)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
c("P30.LPS.MIG","P30.LPS.TDT"),
c("P30.LPS.CTR","P30.LPS.TDT"),
c("P30.PBS.CTR","P30.PBS.MIG"),
c("P30.PBS.MIG","P30.PBS.TDT"),
c("P30.PBS.CTR","P30.PBS.TDT"),
c("P30.LPS.CTR","P30.PBS.CTR")),
label.y = c(0.73,0.78,0.87,0.55,0.65,0.75,0.93),
size=2.35
)
pp.feat.10b2
pp.feat.10c1 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"),
group.by = "new_clusters", #cols = color.cnt,
ncol = 1, pt.size = 0, raster = F) + NoLegend()& geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.feat.10c1
pp.feat.10c2 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"),
group.by = "new_clusters", #cols = color.cnt,
ncol = 1, pt.size = 0, raster = F) + NoLegend() & geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55)&
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55)
pp.feat.10c2
pp.feat.10c3a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)),
features = c("LPS.sig_score"),
group.by = "cluster_cnt", cols = rep(color.cnt[1:3],3),
ncol = 1, pt.size = 0, raster = F)+ NoLegend() & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0,0.7)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = plot.contrast.a,
label.y = (plot.labely.a-1)/10.5,
size=1.75
)
pp.feat.10c3a
pp.feat.10c3b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)),
features = c("LPS.sig_score"),
group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
ncol = 1, pt.size = 0, raster = F) + NoLegend() & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0.35,1)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = plot.contrast.b,
label.y = (plot.labely.b+1)/8.6,
size=1.5
)
pp.feat.10c3b
plot11
DAM score
DAM.sig <- read.csv("../../../20230811_10x_LYN/analysis_plus_exogene/figures1002/new/ranking_of_DAM_indicator_genes.csv")
DAM.sig[1:8,]## 锘縍anking Gene Frequence Total.score
## 1 1 Spp1 11 35.98221
## 2 2 Apoe 11 31.44704
## 3 3 Ifitm3 9 20.82901
## 4 4 Ifi27l2a 9 20.44047
## 5 5 Lyz2 11 20.22463
## 6 6 Cst7 9 19.48372
## 7 7 Cd74 7 18.24720
## 8 8 Lgals3 7 18.24180DAM.list <- list(top50=DAM.sig$Gene[1:50],
top100=DAM.sig$Gene[1:100],
top250=DAM.sig$Gene[1:250],
top500=DAM.sig$Gene[1:500])DAM.list## $top50
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a" "Lyz2" "Cst7"
## [7] "Cd74" "Lgals3" "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1" "H2-Aa" "Fxyd5"
## [19] "Ccl4" "Cybb" "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3" "Axl" "Anxa5"
## [31] "Isg15" "Lgals1" "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5" "Il1b" "Crip1"
## [43] "Slfn2" "B2m" "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1"
##
## $top100
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a" "Lyz2" "Cst7"
## [7] "Cd74" "Lgals3" "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1" "H2-Aa" "Fxyd5"
## [19] "Ccl4" "Cybb" "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3" "Axl" "Anxa5"
## [31] "Isg15" "Lgals1" "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5" "Il1b" "Crip1"
## [43] "Slfn2" "B2m" "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1" "Fcgr4" "Il1rn" "Wfdc17" "Cxcl2"
## [55] "Cxcl16" "Prdx1" "Rtp4" "H2-D1" "Pkm" "Stat1"
## [61] "Anxa2" "Gpnmb" "Zbp1" "Ftl1" "Ldha" "Npc2"
## [67] "Ly6a" "Oasl2" "Gapdh" "Prdx5" "Gbp2" "Grn"
## [73] "Ifi207" "Ifitm2" "Tlr2" "Txn1" "Phf11b" "Ctsz"
## [79] "Pianp" "Cd36" "Itgax" "Fgl2" "Ccl6" "Iigp1"
## [85] "Ccl7" "H2-K1" "Pld3" "Tpi1" "Akr1a1" "Usp18"
## [91] "Rab7b" "Tmsb10" "Cd44" "Aldoa" "Hcar2" "Acod1"
## [97] "Cd300lb" "Ctsc" "Gpr65" "H2-Q7"
##
## $top250
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a" "Lyz2" "Cst7"
## [7] "Cd74" "Lgals3" "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1" "H2-Aa" "Fxyd5"
## [19] "Ccl4" "Cybb" "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3" "Axl" "Anxa5"
## [31] "Isg15" "Lgals1" "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5" "Il1b" "Crip1"
## [43] "Slfn2" "B2m" "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1" "Fcgr4" "Il1rn" "Wfdc17" "Cxcl2"
## [55] "Cxcl16" "Prdx1" "Rtp4" "H2-D1" "Pkm" "Stat1"
## [61] "Anxa2" "Gpnmb" "Zbp1" "Ftl1" "Ldha" "Npc2"
## [67] "Ly6a" "Oasl2" "Gapdh" "Prdx5" "Gbp2" "Grn"
## [73] "Ifi207" "Ifitm2" "Tlr2" "Txn1" "Phf11b" "Ctsz"
## [79] "Pianp" "Cd36" "Itgax" "Fgl2" "Ccl6" "Iigp1"
## [85] "Ccl7" "H2-K1" "Pld3" "Tpi1" "Akr1a1" "Usp18"
## [91] "Rab7b" "Tmsb10" "Cd44" "Aldoa" "Hcar2" "Acod1"
## [97] "Cd300lb" "Ctsc" "Gpr65" "H2-Q7" "Cdkn1a" "Psat1"
## [103] "Trim30a" "Cxcl14" "Srgn" "Mmp12" "Bcl2a1b" "Tap1"
## [109] "AB124611" "Xaf1" "Ly6e" "Psme1" "Ctsl" "Sp100"
## [115] "Cxcr4" "Psmb8" "AA467197" "Pgk1" "Emp3" "Csf1"
## [121] "Mcemp1" "Gusb" "Pim1" "Ctse" "Cox4i1" "Il12b"
## [127] "Msrb1" "Npm1" "Alcam" "Psme2" "Apoc2" "Bhlhe40"
## [133] "Stmn1" "Gm4951" "Pla2g7" "Plaur" "Tor3a" "Hspe1"
## [139] "Tpt1" "Ndufa1" "Flt1" "Tgfbi" "Cox6c" "Irgm1"
## [145] "Ifit2" "Uba52" "Igf2r" "Bola2" "Ank" "Tyrobp"
## [151] "Tpm4" "Ass1" "Ms4a4c" "Ifit1" "Ybx1" "Sod2"
## [157] "Cox8a" "Fam20c" "Oas1a" "Arg1" "Ms4a7" "Smpdl3a"
## [163] "Sh3pxd2b" "Fau" "Gnas" "Phf11d" "Ehd1" "Saa3"
## [169] "Cox5a" "Atox1" "Id3" "Lrpap1" "Olr1" "Sh3bgrl3"
## [175] "Sash1" "Hint1" "Il2rg" "Ctsd" "Postn" "H2-T23"
## [181] "Ucp2" "Sdc3" "Uqcrq" "Cox6a1" "Cd300lf" "Syngr1"
## [187] "Timp2" "Atp5e" "Id2" "S100a6" "Cd14" "Tubb6"
## [193] "Anp32b" "Fcgr2b" "Cd83" "Psmb9" "Bcl2a1a" "Aprt"
## [199] "Mfsd12" "Psap" "Cox6b1" "Lilr4b" "Plac8" "Glrx"
## [205] "Scimp" "Rilpl2" "Psmb6" "Gm11808" "Chmp4b" "Actr3b"
## [211] "Ly86" "Fundc2" "Ifi211" "Hif1a" "Serf2" "Dram1"
## [217] "Parp14" "Ptgs2" "Cxcl9" "Myo5a" "Gabarap" "Cyp4f18"
## [223] "Shisa5" "Lilrb4a" "Cpd" "Iqgap1" "Slfn5" "Tnfaip2"
## [229] "Nme1" "Cotl1" "Tagln2" "Mt1" "Mpeg1" "C3"
## [235] "Ube2l6" "Clic4" "Naaa" "Gas7" "Sdcbp" "Nampt"
## [241] "Ell2" "Samhd1" "Rtcb" "Eef1g" "Hmgb2" "Gng5"
## [247] "Nfil3" "Adora1" "Tpd52" "Ptger4"
##
## $top500
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a"
## [5] "Lyz2" "Cst7" "Cd74" "Lgals3"
## [9] "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1"
## [17] "H2-Aa" "Fxyd5" "Ccl4" "Cybb"
## [21] "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3"
## [29] "Axl" "Anxa5" "Isg15" "Lgals1"
## [33] "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5"
## [41] "Il1b" "Crip1" "Slfn2" "B2m"
## [45] "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1" "Fcgr4" "Il1rn"
## [53] "Wfdc17" "Cxcl2" "Cxcl16" "Prdx1"
## [57] "Rtp4" "H2-D1" "Pkm" "Stat1"
## [61] "Anxa2" "Gpnmb" "Zbp1" "Ftl1"
## [65] "Ldha" "Npc2" "Ly6a" "Oasl2"
## [69] "Gapdh" "Prdx5" "Gbp2" "Grn"
## [73] "Ifi207" "Ifitm2" "Tlr2" "Txn1"
## [77] "Phf11b" "Ctsz" "Pianp" "Cd36"
## [81] "Itgax" "Fgl2" "Ccl6" "Iigp1"
## [85] "Ccl7" "H2-K1" "Pld3" "Tpi1"
## [89] "Akr1a1" "Usp18" "Rab7b" "Tmsb10"
## [93] "Cd44" "Aldoa" "Hcar2" "Acod1"
## [97] "Cd300lb" "Ctsc" "Gpr65" "H2-Q7"
## [101] "Cdkn1a" "Psat1" "Trim30a" "Cxcl14"
## [105] "Srgn" "Mmp12" "Bcl2a1b" "Tap1"
## [109] "AB124611" "Xaf1" "Ly6e" "Psme1"
## [113] "Ctsl" "Sp100" "Cxcr4" "Psmb8"
## [117] "AA467197" "Pgk1" "Emp3" "Csf1"
## [121] "Mcemp1" "Gusb" "Pim1" "Ctse"
## [125] "Cox4i1" "Il12b" "Msrb1" "Npm1"
## [129] "Alcam" "Psme2" "Apoc2" "Bhlhe40"
## [133] "Stmn1" "Gm4951" "Pla2g7" "Plaur"
## [137] "Tor3a" "Hspe1" "Tpt1" "Ndufa1"
## [141] "Flt1" "Tgfbi" "Cox6c" "Irgm1"
## [145] "Ifit2" "Uba52" "Igf2r" "Bola2"
## [149] "Ank" "Tyrobp" "Tpm4" "Ass1"
## [153] "Ms4a4c" "Ifit1" "Ybx1" "Sod2"
## [157] "Cox8a" "Fam20c" "Oas1a" "Arg1"
## [161] "Ms4a7" "Smpdl3a" "Sh3pxd2b" "Fau"
## [165] "Gnas" "Phf11d" "Ehd1" "Saa3"
## [169] "Cox5a" "Atox1" "Id3" "Lrpap1"
## [173] "Olr1" "Sh3bgrl3" "Sash1" "Hint1"
## [177] "Il2rg" "Ctsd" "Postn" "H2-T23"
## [181] "Ucp2" "Sdc3" "Uqcrq" "Cox6a1"
## [185] "Cd300lf" "Syngr1" "Timp2" "Atp5e"
## [189] "Id2" "S100a6" "Cd14" "Tubb6"
## [193] "Anp32b" "Fcgr2b" "Cd83" "Psmb9"
## [197] "Bcl2a1a" "Aprt" "Mfsd12" "Psap"
## [201] "Cox6b1" "Lilr4b" "Plac8" "Glrx"
## [205] "Scimp" "Rilpl2" "Psmb6" "Gm11808"
## [209] "Chmp4b" "Actr3b" "Ly86" "Fundc2"
## [213] "Ifi211" "Hif1a" "Serf2" "Dram1"
## [217] "Parp14" "Ptgs2" "Cxcl9" "Myo5a"
## [221] "Gabarap" "Cyp4f18" "Shisa5" "Lilrb4a"
## [225] "Cpd" "Iqgap1" "Slfn5" "Tnfaip2"
## [229] "Nme1" "Cotl1" "Tagln2" "Mt1"
## [233] "Mpeg1" "C3" "Ube2l6" "Clic4"
## [237] "Naaa" "Gas7" "Sdcbp" "Nampt"
## [241] "Ell2" "Samhd1" "Rtcb" "Eef1g"
## [245] "Hmgb2" "Gng5" "Nfil3" "Adora1"
## [249] "Tpd52" "Ptger4" "Eif2ak2" "Areg"
## [253] "Rsad2" "Slc31a1" "Gm2000" "Cox7c"
## [257] "Tmem256" "Cox5b" "Cyba" "Il18bp"
## [261] "Selenow" "Myl6" "H2-DMa" "Rai14"
## [265] "Dab2" "Hmox1" "Tmsb4x" "Ndufa2"
## [269] "Cd68" "Ccdc86" "Atp6v1a" "Cox7a2"
## [273] "Calm1" "Uqcrh" "Socs2" "Gpi1"
## [277] "Ubl5" "Colec12" "Ifit3b" "Scpep1"
## [281] "S100a11" "F3" "Ms4a6d" "Hacd2"
## [285] "Chil3" "Tuba1c" "Rap2b" "Gp49a"
## [289] "Milr1" "Fcgr1" "Stat2" "Atp5j2"
## [293] "Uqcr10" "Ssr4" "Bcar3" "Gm49368"
## [297] "Tmem106a" "Tubb5" "Naca" "Hspa8"
## [301] "Atp5k" "Bag1" "Sec61b" "Gyg"
## [305] "Cox7b" "Ly6c2" "Top2a" "Aldh2"
## [309] "Dpp7" "Eif3k" "Cd48" "C4b"
## [313] "Pycard" "Atp5l" "Uqcr11" "Osm"
## [317] "Prelid1" "Rnf213" "Prdx4" "Arpc1b"
## [321] "Ndufv3" "Sp110" "Ndufa13" "Abca1"
## [325] "Gm1673" "Wdr89" "Sh3bgrl" "Smdt1"
## [329] "Gatm" "Gpr84" "Slamf8" "Ccrl2"
## [333] "Tomm7" "Gas8" "Ly6i" "Bcl2a1d"
## [337] "Esd" "Dhx58" "Ctsa" "Rxrg"
## [341] "Prdx6" "Ndufc1" "Polr2f" "Sdc4"
## [345] "Atp5j" "Litaf" "Atp6v0e" "Cspg4"
## [349] "Ranbp1" "Ifi35" "Fcer1g" "Calm3"
## [353] "Rhoc" "Pde4b" "Atp5g1" "Gpx3"
## [357] "Psmb1" "Gpx1" "Eef1a1" "Ly9"
## [361] "Igtp" "H2-Q6" "Herc6" "Cd9"
## [365] "Ran" "Hebp1" "Eno1" "Cdk18"
## [369] "Eef1b2" "Gbp7" "Glipr1" "Atp6v1f"
## [373] "H2-DMb1" "Btf3" "Slc25a3" "Brdt"
## [377] "Csf2ra" "Eif3f" "Hpse" "Gm14305"
## [381] "Gm14410" "H2-Q4" "Ndufb9" "Epsti1"
## [385] "Dnaaf3" "Pf4" "S100a4" "Atp5h"
## [389] "Apobec3" "Hsp90ab1" "Tnf" "Ctss"
## [393] "Gas6" "Gbp3" "Gng12" "Nceh1"
## [397] "Mgst1" "Cfl1" "Dtnbp1" "Slc2a6"
## [401] "Eif5a" "Atp5c1" "ptchd1" "Ptchd1"
## [405] "Psma7" "Lap3" "Rbm3" "Cycs"
## [409] "Cox6a2" "Abcg1" "Prr15" "Ahnak"
## [413] "Ndufb7" "Nrp1" "Lmnb1" "Ninj1"
## [417] "Mrpl33" "Arpc3" "Phlda1" "Acp5"
## [421] "Atp5g3" "C5ar1" "Arl5c" "Parp9"
## [425] "Ndufb8" "Txndc17" "Tbca" "Gm49339"
## [429] "Tma7" "Fblim1" "Eif3h" "Psme2b"
## [433] "Mrpl30" "Arpc2" "Ptma" "Rac2"
## [437] "Ctsh" "Dcstamp" "Clec4n" "Dbi"
## [441] "H2-T22" "Sem1" "Msr1" "Bax"
## [445] "S100a10" "Fabp3" "Snrpg" "Bcl2a1"
## [449] "Serpine2" "Snrpd2" "Cd180" "Pgam1"
## [453] "2310061I04Rik" "S100a1" "Serpinb6a" "Actg1"
## [457] "Uqcc2" "Tuba1b" "Neurl2" "Gcnt2"
## [461] "Smc2" "Atp5b" "Ifih1" "Nap1l1"
## [465] "Cd274" "Rgs1" "Parp12" "Serpine1"
## [469] "Nsa2" "Cebpb" "Csf2rb" "Cfb"
## [473] "Ndufa6" "Sdf2l1" "Anxa4" "Psma5"
## [477] "Nfkbiz" "Ctla2a" "Atp5o" "Ube2l3"
## [481] "Lipa" "Pfn1" "Ndufa7" "Ndufs6"
## [485] "Psmb10" "Mapkapk2" "Aif1" "Smc4"
## [489] "Itgb1bp1" "Nme2" "Pfdn5" "Rassf3"
## [493] "1810058I24Rik" "Pgls" "Clec4d" "Mrpl54"
## [497] "Tmem160" "Pomp" "Slc7a11" "F10"lapply(DAM.list, length)## $top50
## [1] 50
##
## $top100
## [1] 100
##
## $top250
## [1] 250
##
## $top500
## [1] 500GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top50, "DAM.sig_top50")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DAM.list$top100, "DAM.sig_top100")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DAM.list$top250, "DAM.sig_top250")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DAM.list$top500, "DAM.sig_top500")## Summarizing datapp.feat.11a1 <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
raster = T, pt.size = 3.5#,keep.scale = "all"
)
pp.feat.11a1
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
pp.feat.11a2 <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
raster = T, pt.size = 3.5,
keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
pp.feat.11a2
pp.feat.11b1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.feat.11b1
pp.feat.11b2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
same.y.lims = T,
group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
c("P30.LPS.MIG","P30.LPS.TDT"),
c("P30.LPS.CTR","P30.LPS.TDT"),
c("P30.PBS.CTR","P30.PBS.MIG"),
c("P30.PBS.MIG","P30.PBS.TDT"),
c("P30.PBS.CTR","P30.PBS.TDT"),
c("P30.LPS.TDT","P30.PBS.TDT")),
label.y = c(0.6,0.75,0.9,0.6,0.75,0.9,1.05),
size=2.35
)
pp.feat.11b2
pp.feat.11c1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
group.by = "new_clusters", #cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.feat.11c1
pp.feat.11c2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
group.by = "new_clusters", #cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55)
pp.feat.11c2
pp.feat.11c3a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)),
features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
group.by = "cluster_cnt", cols = rep(color.cnt[1:3],3),
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(-0.4,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = plot.contrast.a,
label.y = (plot.labely.a-1)/10.5,
size=1.75
)
pp.feat.11c3a
pp.feat.11c3b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)),
features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(-0.25,1.3)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = plot.contrast.b,
label.y = (plot.labely.b-1)/5.5,
size=1.5
)
pp.feat.11c3b
#saveRDS(GEX.seur,"./Spp1tdt.subset2_P30LPS_P30PBS..new1008.rds")
#GEX.seur <- readRDS("./Spp1tdt.subset2_P30LPS_P30PBS..new1008.rds")newlist1009
dotplot and violin in conditions
GEX.seur## An object of class Seurat
## 17900 features across 15536 samples within 1 assay
## Active assay: RNA (17900 features, 1437 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapmarker1009 <- as.vector(unlist(read.table("./dotplot.subset2_newlist1009.txt")))
marker1009 <- unique(marker1009)
marker1009## [1] "Ms4a6c" "Fn1" "Apoe" "Ms4a6d" "Spp1" "Ms4a4c" "Lilrb4a"
## [8] "C4b" "Tspo" "Ms4a6b" "Axl" "Lpl" "Cxcl16" "Ccl2"
## [15] "Ccl12" "Saa3" "Ccl5" "Il1b" "Ccl3" "Cxcl10" "Ccl4"
## [22] "Il1a" "Tnf" "Fcrls" "Gpr34" "Siglech" "P2ry12" "P2ry13"
## [29] "Tmem119" "Sall1" "Slc2a5" "Cst3"length(marker1009)## [1] 32sum(marker1009 %in% markers.new$gene)## [1] 31marker1009.sort <- (markers.new %>% group_by(cluster) %>%
filter(gene %in% marker1009) %>%
ungroup() %>%
#arrange(desc(avg_log2FC*pct.1),gene) %>%
arrange(desc(avg_log2FC),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$genelength(marker1009)## [1] 32length(marker1009.sort)## [1] 31setdiff(marker1009,marker1009.sort)## [1] "Fcrls"marker1009.new <- c(marker1009.sort[1:8],"Fcrls",marker1009.sort[9:31])pp.dot.x1a <- DotPlot(GEX.seur, features = marker1009, group.by = "new_clusters"
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters")
pp.dot.x1a
pp.dot.x1a <- DotPlot(GEX.seur, features = marker1009.new, group.by = "new_clusters",cols = c("lightgrey","red")
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters")
pp.dot.x1a
pp.dot.x1a + scale_color_gradientn(colours = rev(mapal))
pp.dot.x1b <- DotPlot(GEX.seur, features = marker1009.new, group.by = "cnt",cols = c("lightgrey","red")
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Conditions")
pp.dot.x1b
pp.dot.x1b + scale_color_gradientn(colours = rev(mapal))
pp.violin.x1b <- VlnPlot(GEX.seur, features = marker1009.new, group.by = "cnt", cols = color.cnt,
ncol = 5, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.x1b
pp.violin.x1c <- VlnPlot(GEX.seur, features = marker1009.new, group.by = "cnt", cols = color.cnt,
ncol = 5, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0,7.2)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
c("P30.LPS.MIG","P30.LPS.TDT"),
c("P30.LPS.CTR","P30.LPS.TDT"),
c("P30.PBS.CTR","P30.PBS.MIG"),
c("P30.PBS.MIG","P30.PBS.TDT"),
c("P30.PBS.CTR","P30.PBS.TDT"),
c("P30.LPS.TDT","P30.PBS.TDT")),
label.y = c(6.3,6.8,7.3,6.3,6.8,7.3,7.8)-1,
size=1.75
)
pp.violin.x1c
forGEO
GEX.seur <- readRDS("./Spp1tdt.subset2_P30LPS_P30PBS..new1008.rds")
GEX.seur## An object of class Seurat
## 17900 features across 15536 samples within 1 assay
## Active assay: RNA (17900 features, 1437 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapGEX.seur@meta.data[,grep("snn|pANN",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGGTCGACA-1 Spp1tdt 5001 2062 2.699460 12.437512
## AAACCCAAGGTTCTTG-1 Spp1tdt 6501 2457 2.553453 12.305799
## AAACCCACAAACTCGT-1 Spp1tdt 4188 1841 2.960840 8.572111
## AAACCCACACGCGGTT-1 Spp1tdt 4524 1913 1.215738 9.858532
## AAACCCACATCCGTGG-1 Spp1tdt 2337 1225 3.979461 7.659392
## AAACCCAGTATCGAAA-1 Spp1tdt 4318 1955 2.292728 11.741547
## FB.info S.Score G2M.Score Phase seurat_clusters
## AAACCCAAGGTCGACA-1 P30.LPS.CTR 0.0061184939 -0.06523366 S 0
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT1 0.0002076412 -0.07102627 S 0
## AAACCCACAAACTCGT-1 P30.PBS.MIG 0.0192137320 0.06031440 G2M 7
## AAACCCACACGCGGTT-1 P30.LPS.TDT2 -0.0535299003 -0.05464656 G1 0
## AAACCCACATCCGTGG-1 P30.PBS.TDT 0.0489756368 -0.04211671 S 2
## AAACCCAGTATCGAAA-1 P30.LPS.MIG -0.0314645626 0.02183367 G2M 0
## DoubletFinder0.05 DoubletFinder0.1 sort_clusters age
## AAACCCAAGGTCGACA-1 Singlet Singlet 0 P30.LPS
## AAACCCAAGGTTCTTG-1 Singlet Singlet 0 P30.LPS
## AAACCCACAAACTCGT-1 Singlet Singlet 7 P30.PBS
## AAACCCACACGCGGTT-1 Singlet Singlet 0 P30.LPS
## AAACCCACATCCGTGG-1 Singlet Singlet 2 P30.PBS
## AAACCCAGTATCGAAA-1 Singlet Singlet 0 P30.LPS
## cnt preAnno FB.new LPS.TDT.up LPS.TDT.dn
## AAACCCAAGGTCGACA-1 P30.LPS.CTR Microglia P30.LPS.CTR 0.2671683 0.24670044
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT Microglia P30.LPS.TDT1 0.2190005 0.21501295
## AAACCCACAAACTCGT-1 P30.PBS.MIG Microglia P30.PBS.MIG 0.3136800 0.02545367
## AAACCCACACGCGGTT-1 P30.LPS.TDT Microglia P30.LPS.TDT2 0.2974201 0.09338095
## AAACCCACATCCGTGG-1 P30.PBS.TDT Microglia P30.PBS.TDT 0.2758054 -0.03406151
## AAACCCAGTATCGAAA-1 P30.LPS.MIG Microglia P30.LPS.MIG 0.2300058 0.11176488
## new_clusters age1 cluster_cnt DAM.sig_top50
## AAACCCAAGGTCGACA-1 C7 P30.LPS C7_P30.LPS.CTR 0.12647123
## AAACCCAAGGTTCTTG-1 C7 P30.LPS C7_P30.LPS.TDT 0.04136910
## AAACCCACAAACTCGT-1 C2 P30.PBS C2_P30.PBS.MIG -0.03826143
## AAACCCACACGCGGTT-1 C7 P30.LPS C7_P30.LPS.TDT 0.17479143
## AAACCCACATCCGTGG-1 C2 P30.PBS C2_P30.PBS.TDT 0.06229684
## AAACCCAGTATCGAAA-1 C7 P30.LPS C7_P30.LPS.MIG 0.19221428
## DAM.sig_top100 DAM.sig_top250 DAM.sig_top500 LPS.sig_score
## AAACCCAAGGTCGACA-1 0.09610184 0.09546331 0.2497200 0.5988063
## AAACCCAAGGTTCTTG-1 0.02314892 0.06780191 0.2634468 0.6298800
## AAACCCACAAACTCGT-1 -0.02304032 -0.01821938 0.1498645 0.2308747
## AAACCCACACGCGGTT-1 0.26609309 0.19238189 0.3066293 0.6339013
## AAACCCACATCCGTGG-1 -0.06607385 -0.05013960 0.1023302 0.1488521
## AAACCCAGTATCGAAA-1 0.07223828 0.07965503 0.2186862 0.4894010#saveRDS(GEX.seur,"I:/Shared_win/projects/202310_Spp1DAM/forGEO/seur_obj/Spp1TDT_P30LPS.final.seur_obj.rds")